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ABSTRACT 

Monte Carlo approaches to radiation transport have several attractive properties such as simplicity of imple- 
mentation, high accuracy, and good parallel scaling. Moreover, Monte Carlo methods can handle complicated 
geometries and are relatively easy to extend to multiple spatial dimensions, which makes them potentially in- 
teresting in modeling complex multi-dimensional astrophysical phenomena such as core-collapse supernovae. 
The aim of this paper is to explore Monte Carlo methods for modeling neutrino transport in core-collapse su- 
pernovae. We generalize the implicit Monte Carlo photon transport scheme of Fleck & Cummings and gray 
discrete-diffusion scheme of Densmore et al. to energy-, time-, and velocity-dependent neutrino transport. 
Using our ID spherically-symmetric implementation, we show that, similar to the photon transport case, the 
implicit scheme enables significantly larger timesteps compared with explicit time discretization, without sac- 
rificing accuracy, while the discrete-diffusion method leads to significant speed-ups at high optical depth. Our 
results suggest that a combination of spectral, velocity-dependent, implicit Monte Carlo and discrete-diffusion 
Monte Carlo methods represents a robust approach for use in neutrino transport calculations in core-collapse 
supernovae. Our velocity-dependent scheme can easily be adapted to photon transport. 

Subject headings: Hydrodynamics, Neutrinos, Radiative Transfer, Stars: Evolution, Stars: Neutron, Stars: 
Supernovae: General 



1. INTRODUCTION 

Core-collapse supernovae (CCSNe) are among the most en- 
ergetic explosions in the Universe. They mark the end of mas- 
sive star evolution and are powered by the release of gravita- 
tional energy in the collapse of the stellar core to a proto- 
neutron star (PNS). Despite decades of effort, the details 
of the explosion mechanism remain obscure and represent a 
formidable computational challenge. Simulations in spheri- 
cal symmetry with the latest nuclear and neutrino physics, so- 
phisticated neutrino transport, and up-to-date progenitor mod- 
els fail to explode, suggesting that multi-dimen sional effects 
are probably crucial for producing explosions (iHerant et al.l 
[19921 [1991 iBurrows et al.l[T99l Llanka & Muellerill996h . In- 
deed, modern 2D (axisymmetric) simulations, while still am- 
biguous and problematic, exhibit fluid instabilities and tur- 
bulen ce that lead to more f avorable conditions for explo- 
sion ( iMarek & Jankal 120091: lOtt et al.l 120081: lYakunin et al 

m Moreover, re cent ca lculations by [Nordhaus et al. 
: lmiwaki et alJ (l2012h : lHanke et alJ (1201 Ih show that 
the role of the third spatial dimension cannot be neglected, 
and conclusive CCSN simulations will have to be carried out 
in fufl 3D. 
One of the most important ingredients in modeling CCSNe 
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is neutrino transfer Neutrinos play a crucial role in trans- 
porting energy from the PNS to the material behind the su- 
pernova shock, influencing the hydrodynamic and thermody- 
namic conditions of the explosion. At the same time, accurate 
neutrino transport is one of the most complicated and compu- 
tationally expensive aspects of numerical CCSN modeling. 

The transport methods used in previous ID and 2D sim- 
ulations of CCSNe exhibit drawbacks that are likely to be- 
come particularly pronounced in 3D calcul ations. For ex- 
ample, the ray-by-ray me t hod (used, e.g., in j Marek & Jankal 
l2009KBruennetal.ll2006l:lTakiwakietal.ll2012b solves a se- 
ries of coupled ID transport calculations along a number 
of radial rays. While computationally less expensive com- 
pared with a full 3D scheme, this method does not incor- 
porate lateral transport, exaggerates local heating and cool- 
ing, and cannot easily follow off-center motions . The ^Af 
scheme (used, e.g., in 'Liebendo rfer et all l2004t lOtt et al.l 
120081: ISumivoshi &Y amada 2012), while adequate for ID 
calculations, suffer s from so-call ed ray-effects in the higher 
dimensional case dCastoiJ 120041) an d involves a complex 
solution and parallelization scheme (iMcClarren et al.l 120081 : 
ISwestyll2006l) . Although im proving on these dra wbacks is the 
topic of ongoing research (iGodoy & Liul 120121) . it is worth- 
while to explore alternative approaches to neutrino transport. 
One such approach is the Monte Carlo method and the aim 
of this paper is to explore its use in core-collapse supernova 
simulations. 

This paper is organized as follows. In Section |2l we 
summarize the current status of the CCSN simulations, af- 
ter which we present a more detailed introduction to Monte 
Carlo transport methods methods (Section |3). Then, in Sec- 
tion m we describe a simple Monte Carlo method for solv- 
ing the equations of time-dependent radiative transfer For 
this, we restrict ourselves to the ID spherically-symmetric 
problem with a static matter background that is assumed 
to emit, absorb, and scatter radiation In Section |5] we de- 
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scribe some key aspects of a widely used method for the 
time dis cretization of the nonlinear photon transport equa- 
tions by iFleck & Cummingsl (Il971h . In Section |6l we ex- 
tend this method to neutrino transport and provide a Monte 
Carlo interpretation of the resulting equations. In Section |7J 
w e generalize the discr ete-diffusion Monte Carlo scheme 
of lDensmore et al.l(l2007h to energy-dependent neutrino trans- 
port. In Section [8] we describe the extension of this scheme 
to the case when matter is moving. In Section |9] we present 
tests of the numerical implementation of these schemes, while 
in Section [To] we provide conclusions and thoughts about fu- 
ture work. 

2. SUMMARY OF CURRENT CORE-COLLAPSE SUPERNOVA 
SIMULATIONS 

Some basic aspects of the CCSN mechanism are well es- 
tablished. The collapse of the evolved stellar core to a PNS 
and its evolution to a compact cold neutron star (NS) involve 
huge amounts of gravitational energy (~3x lO^-' ergs). The 
explosion mechanism must convert a ~ 10^ '-erg fraction of 
this energy into the kinetic and internal energy of the explod- 
ing stellar envelope to match observations of core-collapse su- 
pemovae. However, after four decades of research, the details 
of this process remain obscure. 

The hydrodynamical shock wave produced by core bounce 
stalls soon after formation, an d it must be reenergized to 
lead to a supernova explosion (iBethd [19901) . The delayed 
neutrino mechanism relies on an imbalance between neu- 
trino heating and cooling behind the shock to deposit suffi- 
cient energy to revive the shock and drive the explosion on a 
timescale of hundreds of milliseconds. However, in spher- 
ical symmetry, this m echanism has been shown to fail for 
regular massive stars ([Burrows et al][1995t iRampp & Janki 
20W; 'Liebendorfer et al."2001l [2005t [Thompson et al.|[2003| 
Sumiyoshi et al. 2005), while fo r the l ow-mass S.SM ^ pro- 
genitor of lNomoto & Hashimotol (Il988h . iKitaura et alj ( l2006h 
obtain a spherical explosion a fter a short post-bounce delay 
(see also lBurrows et al.]|2007ah . This progenitor can explode 
in ID because its envelope is extremely rarefied, but the ex- 
plosion energy is too low (< 10^'^ erg) to match observations 
of typical CCSNe. 

Increases in computer power in the 1990s enabled de- 
tailed numerical simulatio n s in 2 D jBurrows & Frvxelll 
T99l iHerant et all [T99l [T994I: iBurrows et al.l Il995t 
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and potential importance of multi-D hydrodynamical insta- 
bilities and neutrino-driven convection in the core-collapse 
supernova phenomenon. More recent calculations in 2D have 
shown that these instabilities and con vection increase the 
dwell time of matter in the gain region (iMurphv & Burrows[ 
120081) . a region where neutrino heating exceeds neutrino 
cooling. This results in greater neutrino energy deposi- 
tion efficiency behind the shock and, thus, creates more 
favorable c onditio ns for explosion (iMurphv & Burrows 
2008t [Bimo ws & G oshvl[T99l lJankall2001l: iThompson et al 
iPeic] 



20051; iPeicha & ThompsonI 120 12[) . with some of these 
calculations leading to weak delayed neutrino-driven 
explosions (Buras et al. 2006b a; Bruennetal. 2006", 
Mezzacappa et al. 2 007; Bruenn et al. 2009; Marek & Janka 
[2009; Yak unin et al.l[2010l) ! 

However, despite obtaining explosions, these simulations 
pose new questions. First of all, the explosion energies ob- 
tained in these 2D simulations are typically one or two orders 
magnitude smaller than the canonical value. Moreover, these 



exploding models em ploy a soft version o f the n uclear equa- 
tion of state (EOS) bv lLattimer & Swestv[ (1 199 lb with an in- 
compressibility at nuclear densities, K, of 180 MeV. Such a 
soft EOS is r uled out by the recen t observation of a ^2Mc^ 
neutron star (iDemorest et al.l 120101) . iMarek & Jankal (120091) 
did n ot obtain an explosion for their model with K = 263 
MeV dHillebrandt & Wolfj[l985h . suggesting that it may be 
harder to obtain explosions with stiff er EOSs. Furthermore, 
while the Parching and Tokyo groups have f ound marginal 
explosions (iMarek & Jankall2009i: [siiwa et al.ll2010h . the Oak 
Ridge group reports stron ger and earlier explosions for a 
wider range of progenitors ( Bruenn et al. 2009 : Yakuni n et aTl 
[20T0I) . though IBurrows et al.l (12006.. 2007b.) ; . Ott et al.. J2OO8I) 
did not see neutrino-driven explosions for progenitors greater 
in mass than ^8.8 M0. 

Although it is not yet entirely clear why 2D simulations 
by different gr oups produc e different results, the marginality 
of explosion in IMarek & Janka (2009[) and lSuwa et al. (20ly) 
hints at the possible importance of the third spatial dimension 
in explosion dynamics. Three-dimensional (3D) fluid dynam- 
ics has different flow patterns than in 2D. This fact could have 
an impact on the existence and the growth rate of nonradial 
hydrodynamic instabilities in the supernova core, which could 
alter the dynamics of the neutrino-driven ex p losion. Indeed, 
recent simu lations bv [ Nordhaus et al.l (120101) . iTakiwaki et al.l 
( I2OI2I) . and lHanke et alj ( 1201 ll) found significant differences 
in the explosion dynamics between 2D and 3D simulations. 

3. DETERMINISTIC AND MONTE CARLO TRANSPORT 

Two fundamentally different computational approaches ex- 
ist to solve the radiation transport equations, each with well- 
established schools of thought, and with advantages and draw- 
backs. They are the deterministic approach and the Monte 
Carlo approacqj. 

Deterministic methods involve the discretization of the full 
or approximate transport equation on a phase space grid, gen- 
erating a coupled system of algebraic equations. The opti- 
mal way to represent the transport equation on these grids 
for a given situation is frequently far from obvious and is a 
research topic in itself. For example, one may chose finite- 
difference, finite element, or finite volume representations. In 
the momentum-angle variables, discrete ordinates (as in the 
Sn method) or spherical har monic expansi ons (as in the Pn 
method) are often employed ( ICastorl 120041) . Once the equa- 
tions have been discretized, the solution of the resulting sys- 
tem of equations is completely "determined" for given ini- 
tial and boundary conditions. A numerical solution of this 
system produces the global (i.e., over all of phase space) so- 
lution of the transport equation, and provides numerical es- 
timates of the radiation field in the entire problem domain. 
The global nature of such solutions is one of the main ad- 
vantages of deterministic methods. However, the discretiza- 
tion process introduces (often significant) truncation errors; 
for e xample, a simple Pn may suffer from negative ener- 
gies dMcClarren et al.ll2008l). Reducing such error s is an area 
of active research (e.g.. lMcClarren & Hauck|[2010b 

There are various deterministic approac hes to both approx - 
imate [e.g., the diffusion approximation (lPomraninglll973l) 1 
and full multi-angle and multi-energy transport. For the lat- 
ter, one of the most widely used approaches is the discrete- 
ordinates Sn method which solves the transport equation 

' Recently, hybrid methods that combine Monte Carlo and deterministic 
methods have also found some success ( Wollaber & Larsen 20091) . 
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along several particular directions in each spatial zone (ICastoiJ 
120041) . However, such methods have several drawbacks. 
First and foremost, they suffer from ray effects (iMorelet alj 
l2003h . Because of the discrete nature of the angular repre- 
sentation, this method introduces large spatial oscillations in, 
e.g., energy density. Also, Sn methods employ a very com- 
plex solution and parallelization procedures. For large sys- 
tems, direct inversion of the transport operator can be very 
inefficient , forcing one to resort to complicated iterative ap- 
proaches (I Adams & LarsenI 120021) . A further limitation of 
such methods that has emerged more recently i s their some- 
what limited parallel scalability. iSwestvl (120061) shows that a 
variant of the Sn scheme based on Krylov iterative frameork 
can scale well up to only ^ 256 cores for a uniform grid (al- 
though the exact details somewhat depend on problem param- 
eters). This is significant because any 3D radiation transport 
calculation is likely to require parallel calculations on many 
thousands of processors. Improving the para llel scalability 
of such methods is an area of active research (ISwestvl 120061 
iGodov & Liull20T2i) . 

In contrast with deterministic approaches, in Monte Carlo 
methods one does not solve the transport equation; instead, 
such methods employ pseudo-random number sequences to 
directly simulate the transport of radiation particles through 
matter. Since the maximum number of particles that one can 
simulate is constrained by computer memory and CPU power, 
typically many fewer Monte Carlo particles than actual physi- 
cal particles participate in the numerical transport process, im- 
plying that each Monte Carlo particle represents some packet 
of many physical particles. Based on local emissivity, Monte 
Carlo particles are sampled in various zones with various fre- 
quencies, directions, and spatial coordinates. Then, each par- 
ticle is tracked through the problem domain until it crosses 
a boundary or is absorbed. If a sufficiently large number of 
Monte Carlo particles is simulated, then one can obtain an ac- 
curate estimate of the average behavior of the system. This is 
the basic idea behind Monte Carlo radiation transport^- 

In this paper, the term "particle" refers to a single radiation 
particle, such as a photon or a neutrino. We use the term "MC 
particle" or "MCP" to refer to a Monte Carlo particle that rep- 
resents a packet of physical particles. The number of physical 
particles represented by a given MCP will be referred to as 
the weight of the particle. 

Monte Carlo methods have some interesting features that 
may be particularly advantageous in multi-dimensional trans- 
port simulations. First, such methods are generally eas- 
ily adapted to work with complicated geometries, meshes, 
and multiple spatial dimensions. This is because the most 
geometry-dependent aspect of Monte Carlo methods consists 
of the algorithms that deal with tracking MC particles through 
spatial zones in the problem domain. Most of the rest of the 
algorithm represents geometry-independent operations such 
as equation of state and opacity calculations, and calculations 
that depend on the length of particle paths. Once this track- 
ing is implemented for a given mesh type, th e rest of a Mon te 
Carlo algorithm is relatively straightforward (lGentilelt2009l) . 

Second, Monte Carlo methods model physical processes in 
a more direct and simple way than deterministic methods. 
For example, anisotropic scattering is handled easily, just by 

'" IZinkI j2008f) suggested a Monte Carlo discretization of the general 
relativistic transport equations. In this scheme, contrary to traditional 
Monte Carlo radiation transpoit, one solves the transport equations using 
Monte Carlo methods, instead of directly simulating radiation transport using 
pseudo-random numbers. 



changing the Monte Carlo partic le's direction when a scat- 
tering event occurs (ICastoiJ 12004b . Since scattering is mod- 
eled by deflecting the paths of MC particles in a simulation, it 
can represent the angular beh avior of a scattering kernel with 
more fidelity (iGentild 120091) . The angle of the scattered ra- 
diation particle can be chosen from a probability distribution 
function that can easily be constructed for scattering kernels 
of any (physically reasonable) functional form. The same is 
also true for the operation of selecting the energy of scattered 
particles in the case of inelastic scattering. Incorporation of 
anisotropic and inelastic scattering in deterministic methods is 
far more involve d and much less straightforwa rd [but is nev- 
ertheless doable (iMezzacappa & Bruennlll993l) 1. 

The MC method can also be modified to account for 
material motion in a relatively straightforward manner 
using a mixed-frame fo rmalism (iMihalas & Klein] 119821 : 
iHubenv & Burrowsl l20071 Emission takes places in the fluid 
frame and radiation particles are Lorentz-transformed into 
the Eulerian lab frame, where transport is performed. This 
method produces the correct distribution for a fluid moving 
with relativistic velocity. Velocity-dependence in determinis- 
tic methods is again much more complicated to implement. 

Monte Carlo methods also have the advantage that if the 
entire problem domain (i.e., meshes, hydrodynamic and ther- 
modynamic variables, etc.) can fit into the memory of one 
CPU node, then parallelization is trivial and strongly scalable. 
One just simulates copies of the problem on a number of pro- 
cessors, where each processor carries a fraction of the total 
number of particles. Quantities accumulated over all particles 
(e.g., the total emitted or deposited energy or lepton number, 
etc.) are then summed over all proces sors. This app roach is 
usually referred to as mesh replication (lGentilel2009b . Even if 
the problem domain does not fit into the memory of one CPU 
node, it is frequently possible to decompose the domain of the 
problem onto separ ate nodes and maintain a high degree of 
parall el scalability dBrunner et al.ll20()6t iBrunner & BrantlevI 
l2009l) . 

There are, however, also some negative aspects of Monte 
Carlo methods. The most serious property is the noise in- 
trinsic to random processes. Monte Carlo methods exhibit 
statistical fluctuations in quantities such as radiation energy 
density and temperature. According to the central limit theo- 
rem, this statistical error scales as N~^l^, w here A^ is the num- 
ber of MC particles used in the calculation (iKalos & Whitlocfl 
120081) . Because the noise (more rigorously, the standard de- 
viation of calculated quantities) decreases so slowly with the 
number of MC particles, it can take many particles to produce 
a sufficiently smooth solution, and this can make large sim- 
ulations computationally very expensive. Therefore, Monte 
Carlo methods are likely to be more expensive than determin- 
istic methods in lower-dimensional cases. On the other hand, 
Monte Carlo methods are likely to become more efficient with 
increasing number of dimensions, at least for calculations of 
definite int egrals (see, e.g., the discussion in Section IIB of 
IZinkll20()8b . However, numerical integration is not the same 
as direct simulation of radiation transport. Hence, such an ar- 
gument should be taken with a grain of salt. Therefore, it is 
a priori unclear how a Monte Carlo method with good paral- 
lel scaling may fare against deterministic methods in terms of 
total computational cost in multi-dimensional calculations. 

When emission and absorption of radiation leads to non- 
negligible cooling or heating of matter through which radia- 
tion is propagating, then the transport problem becomes «o«- 
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linear. Such a scenario is described by a system of non-linear 
equations with a number of unknowns: the radiation inten- 
sity, the material temperature (here and hereafter we assume 
that material is well-described by temperature, i.e., that it is 
in thermal equilibrium), and the leptonic composition (if we 
are dealing with the transfer of neutrinos with lepton number). 
These equations are coupled due to absorption and emission 
terms - the material cools through emission and heats through 
absorption. Similarly, inelastic scattering also leads to nonlin- 
ear coupling between the material temperature and radiation. 

The state-of-the-art in both deterministic and Monte Carlo 
methods for solving non-linear radiation transport problems 
involves linearizing equations over a timestep and solving 
the resulting linear system during the timestep. Perform- 
ing this linearization produces a linearization error during the 
timestep, but it enables the use of a large portion of the ex- 
isting arsenal of linear transfer methods. Moreover, the lin- 
earization errors can b e mitigated by perfo rming iterations 
within a timestep (e.g.. lBurrows et alll2000l) . Non-linear (or, 
more prec isely, semi-hnea r) solution schemes have also been 
proposed (lN'Kaouall991h . but they are not widely explored in 
practical problems. In this paper, we consider only methods 
that involve a linearization procedure within each timestep. 

The classic and widely used method for nonlinear Monte 
Carlo photon transport is the method of Fleck and Cummings 
(1971, hereafter FC71). This method is known as Implicit 
Monte Carlo (IMC). This method reformulates the nonlinear 
transport equation so that the emission term is treated semi- 
implicitly. This process leads to the effective reduction of 
the emission and absorption opacity, and the appearance of 
a scattering term that effectively replaces a fraction of the ab- 
sorption and re-emission of radiation within a timestep. This 
reduces the coupling between radiation and matter within a 
timestep, enabling much larger timesteps and significantly im- 
proving the stability of the system (Wollaber 2008). Since its 
first publication, the IMC method has successfully been used 
in photon transport, in part becau se of its simplicity, versatil- 
ity, and robustness (lGentilell2009l) . In this paper, we general- 
ize the IMC method to neutrino transport. 

One drawback of the IMC method is that it becomes com- 
putationally inefficient at high optical depth. This is because 
in such regimes the radiation mean-free-path due to effec- 
tive scattering becomes very small, i.e., most of the compu- 
tation is spent in modeling these scatterings. Several methods 
have been suggested to overcome this inefficiency. One of 
the simplest and most efficient such methods is the discrete- 
diffusion Monte Carlo (DDMC) scheme of Densmore et al. 
(2007), developed for the case of gray transport for non- 
moving m atter In this pa p er, we extend the gray DDMC 
scheme of iDensmore et al.l (l2007h to the multi-group case, 
and generalize for moving matter We demonstrate that the 
combination of the IMC scheme at low optical depths with 
the DDMC scheme at high optical depths is an attractive ap- 
proach for neutrino transport in CCSN simulations. 

We stress that in the present work our focus is on Monte 
Carlo neutrino transport, as well as on energy and lepton 
number coupling between radiation and matter The is- 
sue of momentum coupling between radiation and matter is 
not discussed, but is straightforward. The full radiation- 
hydrodynamics scheme and associated simulations will be 
presented in a subsequent publication. 

We point out that neutrino transport in previous time- 
dependent simulations of CCNSe has been performed us- 
ing only deterministic methods. Monte Carlo methods were 



used for the study of neu trino equilibration in static uni- 
form matter (iTubbsl [T978h and for calculations of station- 
ary neutrino transfer in static sph erically symmetric super- 
nova matter (iJanka & HillebraridB [l989l The latter code has 
also been ap plied to the study of neutrino spectrum forma- 
tion (iJanka & Hillebrandt 1989i: iKeil et al. 2003) . neutrino- 
antineutrino annihilation (iJankal Il991h . and for assessing 
the quality of deter ministic transport solvers dJankal [r992l : 
lYamada et al.lll999l) . Unlike these codes, our Monte Carlo 
code is fully time-dependent and can handle energy and lepton 
number coupling between matter and radiation, matter mo- 
tion, as well as diffusion at high optical depth. 

Unless otherwise noted, in the following we use spherical 
polar coordinates and CGS units. 

4. A SIMPLE MONTE CARLO METHOD FOR RADIATION 
TRANSPORT 

In this section, for completeness we describe some salient 
aspects of a simple time-explicit Monte Carlo method for non- 
linear time-dependent radiative transfer For simplicity of il- 
lustration, we consider static matter that emits, absorbs, and 
scatters radiation. Our description closely follows the presen- 
tation in Chapter 3 of Wollaber (2008). We start by writing 
the m ulti-D transport equation for such a system (iPomraningI 
rT97l : 

Idl 

-—ir,n,e,t) + n-'^I{r,n,e,t) = 
c at 

Ka(e,r)[B(£,r)-/(r,ii,£,f)]-Ki(e,r)/(r,n,e,f) 



+ / / >c,{e\n' ^e,n)I(j,n\e',t)dn'de' , (1) 
which is coupled to the material energy equatioio 
p-^{rj) = j j K,{eJ)[I{r,n,e,t)-B{eJ)]dnde 
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--KAe',n' -^ £,n)/(r,n',e',f) 

e' 



->cs{e,n-^ e',n')I{r,n,e,t) 



dildil'dede' , 



(2) 

where / is the radiation specific intensity, t is the time, T is 
the temperature, p is the matter density, c is the speed of light, 
K„ is the total absorption opacity, and k, is the total scattering 
opacity. Also, r is the spatial coordinate, n is a unit vector in 
the radiation particle propagation direction, il is the solid an- 
gle, and M:s{e',n' — > e,n) is the differential scattering opacity 
for scattering from energy and propagation direction {e',n'} 
to {£,n} (for brevity the argument T of function Xs is sup- 
pressed). B is the Planck function if the radiation particles are 
photons (or the Fermi-Dirac function if we are dealing with 
fermions), e is the energy of a single physical radiation parti- 
cle, and Ui„ is the specific internal energy of matter Here and 
hereafter, we define opacity as the inverse mean-free-path of 
radiation particles. 

' ' In the case of neutrinos with lepton number, the transport equation (T) 
is also coupled to the equation for the electron fraction Ye of the material. 
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From here on, we assume that the matter is distributed 
spherically symmetrically and use the spherical polar coor- 
dinate system. Therefore, our system is described by only 
radius {0 < r < R) and one angular variable /i = cos 6, where 
9 is the angle between r and the particle propagation direc- 
tion. The radiative transfer equation in spherical coordinates 
is given by: 



1 9/(r,/i,e,f) dl(r,fi,e,t) l-/i- 9/(r,/i,e,f) 
c dt dr r 



d^i 



= Kaie, T) [B(e, T)-I{x, fi,e,t)] - K,(e, T)I{r, ^, e, t) 

/+1 poo 
/ >ts{£' ,^jl' ~^ £, ^i)I{x,^.' ,£' ,t)d^'de' , (3) 

while the material energy equation (|2]i has the following 
form: 

p—^ix,T) = 2-n- / Ka(e,T)[l(x,fj,,e,t) 

-B{e,T)]dfide + S, (4) 

where function S represents the amount of energy exchange 
between radiation and matter due to inelastic scattering: 



S = (IttY 



OO /-OO fl />1 



JO 



1 J-1 



->c,ie',fj.'^e,fj,)I{x,fj.',e',t) 



-Ks(e,ii^ e',fi')I(x,fi,e,t) 



dede'dfid^' . (5) 



Here again, K^ie' .jj! — > e,/i) is the differential scattering 
opacity for scattering from energy and angle {e', /i'} to {e, /^}. 
The initial conditions are: 



/(r,pi,e,0) = /,(r,/i,e), 
T{r,0) = Ti{r), 



(6) 

(7) 



and the boundary conditions are 

IiR,H,e,0) = M^i,s,t), -1<^<0. (8) 

We assume that our spatial domain r G iQ,R] is split into 
many non-overlapping spatial zones with coordinates r e 
[0-i/2> 0+1/2]' where ; = {!,... ,A^J, ri/2 = 0, and rN,.+i/2 = R- 
The quantities that represent the properties of the matter (such 
as temperature, opacity, emissivity, etc.) are represented on 
these cells within each timestep t„ < t < f„+i with their cell- 
averaged values at f = f„. In the following, we describe a sim- 
ple Monte Carlo method for solving equations (|3]|8]l. 

We start by considering the possible sources and sinks of 
radiation particles that enter the transport equation (|3). For 
instance, in the first timestep, MCPs may be present initially, 
or are born due to the boundary conditions or emission by the 
matter. The energies of the emitted particles are subtracted 
from the material internal energy. By the end of a timestep, 
some MCPs may have been absorbed in the material. The 
energies of these MCPs are added to the material internal en- 
ergy and these MCPs are removed from computer memory. A 
fraction of MCPs may leave the system via the outer bound- 
ary. Other MCPs may continue to exist - these MCPs are 
usually stored in a census in computer memory in preparation 
for the next timestep. At the beginning of the next timestep, 
these MCPs emerge from the census (similar to the situation 



in the first timestep), while boundary conditions and emission 
may supply additional MCPs. Using this synopsis as a guide, 
a natural algorithm emerges with which to perform the Monte 
Carlo sampling procedure: We use random numbers to choose 
the positions, the propagation directions, and the energies of 
the newly-born MCPs. Once this is done, random numbers 
are used to simulate the propagation of these MCPs through 
matter within a timestep. This procedure is described in more 
detail in the following. 

We first choose the weight of MCPs, i.e., the total number 
of radiation particles contained in each MCP For simplicity, 
we assume that each MCP represents A^o radiation particles. 
Each radiation particle within a given MCP has the same po- 
sition, propagation angle, and energy. 

The total number of radiation particles emitted by matter is : 



A/V = Stt^ 



Jo 



-^(eJMeJ) ^,^^^^^^^ ^^^ 



Since each MCP contains A'o radiation particles, the total 
number of MCPs emitted in this process is 



A^7- = RInt(7Vr/A^o 



(10) 



Here Rlnt(jc) is an operator that returns the largest integer that 
is not greater than x plus the quantity K, which is chosen ran- 
domly to be 1 with probability p = {Mr /No}, where the latter 
is the fractional part of Mt/Nq. Otherwise, K is selected to 
be 0. In practice, this is done by sampling a (pseudo) random 
number ^ with uniform distribution on the interval [0,1]. If 
£, < p then K = I, otherwise K = 0. 

The particle's energy in each MCP is chosen according to 
the functional form of Ka{e, T)B{e, T). Due to the isotropy of 
the emission, the MCP angle is chosen uniformly on a unit 
sphere. In practice, this is usually done by choosing p uni- 
formly on the interval [-1,1] from 



p = 2^-l. 



(11) 



Since we use time-centered values of the emissivities within 
any timestep t„ <t < t„+\, the particles are emitted with uni- 
form probability within t S [f„,f„+i]. Hence, the emission time 
of the MCP is chosen as 



t - tn + (tn+l-t„)^. 



(12) 



In order to choose the MCP spatial location, one first recalls 
that for transport problems space is represented by many con- 
nected, non-overlapping spatial zones. We first choose the 
zone in which an MCP is born, after which we select the spa- 
tial location of the MCP within that cell. More specifically, if 
Ntj is the total number of particles emitted in zone j, then an 
MCP is born in that cell with probability Ntj/Nt- Note that 
here we assume that the weight of MCPs is given in terms of 
the total number of physical particles represented by a single 
MCP. If the weight were given in terms of the total energy of 
MCPs, then we would need to use the ratio of the total energy 
of particles emitted in each zone to that of particles emitted in 
all of the zones. For a ID zone defined by [?")-i/2, 0+1/2]' the 
particle location r is chosen according to 



0-1/2 +(0+1/2- 0-1/2) ^C 



1/3 



(13) 



which guarantees a uniform sampling within the cell vol- 
um^. 
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However, in highly-diffusive regimes the sampling of the location of an 
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The number of particles that appear during f G [f„,f„+i] due 
to a boundary source at r = R is obtained by integrating the 
boundary condition ^ over the timestep f G [t„,t„+i], the 
boundary surface areaATiR^, and the angle fi G [-1,0); 



given by 



Nb = RInt 



No 



^ilR(_^i,e,t) 



dtded^b 



(14) 

where the particle location on the boundary, the direction, the 
energy, and its time of emission are selected according to the 
functional forms of Ir. 

During the first timestep, particles may also be present due 
to initial conditions, i.e.: 



Nic = RInt 



cA^o 70 J-i JO 



Ii{r,^,e)r drd ij.de 



(15) 



Here, the particle's cell, the spatial location, direction, and 
energy are again selected randomly, this time using the func- 
tional form of/,. 

Thus, the total number of MCPs contained in the problem 
during the first timestep is 



Ntot =Nt+Nb+ Nic , 
while in subsequent timesteps this number is 

Ntot =Nt+Nb+Nc, 



(16) 



(17) 



where Nc is the number of particles in the census from previ- 
ous timesteps. 

Once an MC particle is introduced into the problem, the 
next task is to transport it through the system (and update all 
the relevant quantities along the way). There are essentially 
three types of events that must be considered and which can 
affect the transport of the MCP: 

• The MC particle could collide with the matter (e.g., 
with an atom, a nucleus, or an electron), 



• the MC particle could leave one cell and enter an adja- 
cent cell with different opacities, or 

• the MC particle could travel without collisions inside 
the cell until the end of the timestep (i.e., while t < t„+i). 

There are three different distances associated with these three 
possibilities: the distance to collision d^, the distance to the 
cell boundary d^, and the distance d, that the particle would 
travel until t = f„+i[3 The distance to the cell boundary can be 
calculated using elementary geometric considerations and is 

MC particle should reflect the gradient of the thermal emissivity within the 
zone, i.e.. particles should be bom with higher probability at po ints within 
the cell where the emissivity is higher jFleck & Canfieldlll984l) . Uniform 
sampling may lead to unphysical results (Densmore 2011, private communi- 
cation). 

'^ In an alternative Monte Carlo approach, the so-called continuous ab- 
sorption method can be used (e.g., FC71). This is explained in more detail in 
Sectionim 



db= < 



'-].1/2-'-'(1-m') 



'-}_1/2-'-'(1-m') 



1/2 



-rfi 



1/2 



-l-r/i 



if ;■ = 1 or 

fj,>0, sm9> ^ 



ifAi<0, sin6i<^ 



(18) 
The distance traveled until the end of the timestep, dt, is sim- 
ply given by 

d, = c(t„^i-t). (19) 

In the simplest case, the distance to a collision can be calcu- 
lated probabiUstically and is given by 



dr =- 






(20) 



where ^ is a random numb er with a unifor m distribution on 
the interval (0, 1] (See, e.g.. lWollabe^l2008l for the derivation 
of this formula.). Once the three distances are calculated, the 
next step is to determine which one is the smallest of the three. 
Depending on which is smallest, the MCP is then moved to 
either the collision location, the cell spatial boundary, or the 
time boundary. Accordingly, the MCP location and time are 
updated using the operation: 



r -^ \/ r- — Irdji + d^ , 



t-^t + d/c, 



(21) 



(22) 



where d = mm{dc,db,dt} is the minimum distance. If d = db, 
then we check whether this boundary is the outer boundary of 
the computational domain. If that is the case, then the MCP 
leaves the system (and, thus, information about the MCP is 
erased from computer memory). Otherwise, the transport 
sampling process begins again in the new spatial zone (with a 
new opacity). If d = d,, the MCP is stored in computer mem- 
ory for the next timestep. 

If d = dc, the type of collision event must be determined. 
The MCP is absorbed with probability of /?„ = Ka/{Ka+Ks) and 
scattered with probability of p., = \-pa = Hs/{na + Ks). If it is 
absorbed, then the MCP energy is deposited into the cell and 
information about the MCP is erased from computer memory. 
If it is scattered, then, once the particle location and time are 
updated according to equations (12111221 1. the new angle and (if 
the scattering is inelastic) energy of the particle are selected 
randomly from the functional form of the scattering kernel. 

Finally, at the end of the timestep, the material temperature 
in each spatial zone is updated according to equation (JDi using 
the information about how many particles (of which energy) 
are emitted, absorbed, or scattered in each zone. This process 
is then repeated for each new timestep, for each MCP. 

4. 1 . The continuous absorption method 

The continuous absorption method is a variance reduction 
mechanism t hat is typically used in practical implementa- 
tions of IMC ( IWollabedl2008h . In this method, one calculates 
four different distances (instead of the three distances in the 
method described above): the distance to the boundary, db, the 
distance traveled by the MCP until the end of the timestep, dt, 
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the distance to scattering, ds, and the distance to absorption, 
da- The distances dj, and dt are again calculated using equa- 
tions ( fTSl l and ( fT9] ), respectively. The distance to scattering is 
calculated probabilistically (similarly to equation l20t: 






(23) 



where ^ is a random number with a uniform distribution on the 
interval (0, 1]. On the other hand, the distance to absorption 
is calculated deterministically in the following way. When an 
MCP propagates a distance dx through a material with absorp- 
tion opacity «;„, then the number of radiation particles N{t) in 
this MCP at time t decreases according to the law 



A^(f)=A^(0)e- 



,r,dx 



(24) 



where A^(0) is the initial number of radiation particles in the 
MCP. An MCP is assumed to be absorbed when only a small 
user-defined fraction <r of the initial radiation particles remains 
in the MCP. The parameter <, is usually chosen to be 0.01 
(FC71). 

The transport algorithm in the continuous absorption 
method is again based on the calculation of the smallest of the 
distances. However, as mentioned above, in this case, we are 
dealing with four different distances: da, d/,, d,,, and d, . If da is 
the smallest of the four, then we deposit all the particle energy 
(and lepton number, if we are dealing with neutrinos with lep- 
ton number) into its spatial cell. If the minimum distance is 
ds, dt, or dt, then we move the MCP to its new location ac- 
cording to equations (ISTT i and (l22l i. We then calculate what 
fraction of the MCP is absorbed according to equation (l24t 
during its propagation to its new location, and deposit the en- 
ergy (and the lepton number, if the particles are neutrinos with 
lepton number) of the absorbed fraction into the spatial cell. 
After that we perform a scattering if d,, is the smallest of the 
four, or move to a new cell if dh is the smallest. 

5. THE FLECK & CUMMINGS METHOD FOR IMPLICIT MONTE 
CARLO PHOTON TRANSPORT 

The first step in almost all of the commonly-used methods 
for solving non-linear transport equations is to linearize them 
over a timestep r„ < f < r„+i. As mentioned in Section[T] this 
linearization introduces discretization errors (that grow with 
the size of the timestep), but it allows use of the large number 
of techniques developed for solving linear radiation transport. 
One of the most well-known and widely-used linearization 
techniques is the Implicit Monte Carlo method suggested by 
FC71. 

Consider a ID spherically symmetric problem with static 
matter that can emit, absorb, and scatter radiation (the gen- 
eralization to multi-D is conceptually trivial). The transport 
equation for such a system is given by equation Q, while the 
material energy equation is given by equation (|4). The IMC 
method reformulates the transport equation (|3) using the ma- 
terial energy equation (|4), so that the emissivity in the former 
equation is treated implicitly. This leads to the appearance of 
two new terms in the transport equation, which look like sink 
and source terms due to some scattering process. This scatter- 
ing is called effective scattering by FC71, and it models ab- 
sorption and re-emission of a photon within a timestep. The 
introduction of effective scattering reduces the stiffness of the 
non-linear coupling between the matter temperature and the 
radiation, significantly improving the stability of the system 
of the equations relative to the case when there is no effective 
scattering (.Larsen & Mercieri.1987,) . 



The central point in this reformulation of the transport equa- 
tion is to approximate the radiation source term KaB using the 
value of the intensity at the current time f e [fn,f,,+i] and using 
the values of the other quantities at the beginning of timestep 
f = f„. In that case, the coupling between the two equations 
would simplify. Specifically, equation ^ can be solved in- 
dependently of equation (|4|i within a timestep, while the re- 
sult of solving equation (O can then be used to solve equa- 
tion (HI within the same timestep. This approximation allows 
for much larger timesteps than the mean absorption and re- 
emission timescale (a very short interval in highly-diffusive 
regions), without compromising accuracjo Although the 
emission term is treated semi-implicitly in the IMC method, 
the term "implicit" is, strictly speaking, a misnomer since the 
rest of the problem parameters must be (explicitly) evaluated 
prior to performing the timestep. Since the original work 
by Fleck & Cummings, the IMC method has been widely 
and succe ssfully used for solving many radiativ e transfer 
problems (lGentilel l200Ill2009; Mc Clarren & Urba tsch 20091 
iKasen et al.ll20IIh . As a prelude to the extension of the IMC 
method to neutrinos in the next section, here we describe 
some of the key aspects of the Fleck & Cummings method 
for photons. 

We start by introducing a new set of variables: 



47r f°° 
Ur=— / Bde, 
c Jo 

4nJ-Bds^ 
_ ir i^aBde 



iTBde ' 



and 



P-- 



1 dUr 

pdU,„ 



(25) 
(26) 

(27) 
(28) 



where Ur is the radiation energy density if in thermodynamic 
equilibrium, and Kp is the Planck mean opacity. Using these 
new functions, we rewrite the transport equation (|3) in the 
following wajH: 

ldl(e,^) dl{e,^) I - ^i^ dl(e , fi) 
c dt dr r dfi 

= KabcUr - KaI(£,lS) - Ks/(e, jS) 

/+1 />oo 
/ Ksie\n' ^e, n)I(p\e')dn'de' , (29) 

while the material energy equation can be transformed into 
the following: 

1 dUr /■' f°° 

— -^ \-KpCUr = 2TT I I Kaldfide + S. (30) 

In most applications, the information (such as tempera- 
ture) is known at the beginning of timestep t = t„ from the 

'"* However, too large t imesteps may lead t o unphysical solu- 
tions jLarsen & Mercieil[l987t IMartin & BrownlEOOTt IDensmore & LarsenI 
[2001 

'^ Here and hereafter, we do not consider changes in Ur, p, and other 
thermodynamic variables due to the motion of matter This issue will be ad- 
dressed in a future publication on full radiation-hydrodynamics calculations. 
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previous timestep or from the initial conditions, and one 
needs to find the solution at the end of timestep t = t„+i. 
We approximate the functions {Ka,Kp,Ks,M:s,b,/3} with con- 
stants {Ka,kp,ks,ks,b,P} that are time-centered values of 
{Ka,Kp,b,(3} within [t„,t„+i]. Obviously, such an approxima- 
tion loses its validity if these functions change rapidly within 
a timestep. In many practical applications, these functions are 
usually given by their values at the beginning of the timestep, 
but, if need be, these can also be extrapolated from their val- 
ues at the previous timestep ("FC71r^ Using this approxima- 
tion in equation ( l30t , we obtain a linear equation for U,'. 

1 dUr [^ f°° 

y;-— + kacUr = 2TT I I kJdfide + S. (31) 

Using this equation, FC7 1 derive an approximate equation for 

t/,.(f)forf e[f„,f„+i]: 



Ur(t) = fnU;, + 2TT 



l-fn 



CK„ 



kal{t)d^de , 



(32) 



1 Jo 



where [/*,, = Urj,+pi^t„S, Af„ = f„+i -f„, S is the time-averaged 
value of S within interval t £ [t„ , t„+i ] : 



S = 



1 
Al, 



S(t)dt, 



and /„ is a new variable defined as 

1 



fn = 



l+aAt„/3ckp 



(33) 



(34) 



where a is a user-defined constant such that a € [0.5, 1] for 
stability. The variable /„ is called the Fleck factor. 

Equation ( |32] | is an important result because Ur{t) at a given 
time t G [f„,f„+i] depends explicitly on the [/,.„ at f = t„ and the 
intensity at f e [f„,f„+i]. The term S can be treated by using 
S from the previous timestep. Hence, if we approximate the 
transport equation (|29^ using {na^Hs^Ki^P} = {ka,ks,>Cs,$} 
and substitute U,- in the resulting equation with the RHS of 
( |32] ). we obtain a transport equation that can be solved inde- 
pendently of the material energy equation (|D within timestep 
t,i <t < f„+i. The resulting transport equation has the follow- 
ing form: 

1 5/(/i, e) dlifi', e') 1 - fi^ dl(fi, e) 

W. — "'"^' — JT — + = 

c at or r 



„bcU: 



dfi 



Jin, e) - Ke/(M, e) - kJin, e) 



+27r 



kab 



kesl(n',s')d^j,'d£' 



1 Jo 



+2tt 



kAs', 



e, /i)/(/i',£')c//i'iie', (35) 



'^ FC71 hint that one can also use time extrapolation to determine temper- 
ature from the values at previous timesteps. However, experience has shown 
that these temperature extrapolations can affect the stability and accuracy of 
the result, especially if the solution method is subject to errors (such as sta- 
tistical noise in an MC calculation). Hence, temperature extrapolation is usu- 
ally avoided in practice, and the problem data are frozen at the beginning of 
timestep iWollabei 2008). Alternatively, one can estimate the temperature at 
the end of the timestep using an additional relatively inexpensive determinis- 
tic calculation ( Wollaber & Larsen 20091) . 



where we have introduced two new variables: 

l^ea — JnKa ■: (3o) 

kes = (i-fn)ka, (37) 

the sum of which equals the total absorption opacity. 

We now explore the physical meaning of the new terms on 
the RHS of equation ( [35] ). Terms keabcU*^^ and keJ look like 
source and sink terms due to emission and absorption of par- 
ticles with absorption opacity n^ [compare these terms to the 
1st and 2nd terms on the RHS of equation (|29l)1. Moreover, 

the terms kesi and 27r-^ Jj /g°° k^Jd^'de' look like sink and 
source terms for scattering. Hence, equation (l35t appears to 
describe the transport of radiation through matter with absorp- 
tion opacity, k^a, and an additional scattering opacity, k^s (in 
addition to k^ and ks). For that reason, in this formalism, a 
portion of true absorption and re-emission within a timestep 
is modeled as an effective scattering process. Parameters kea 
and kes are called by FC71 "effective" absorption and scatter- 
ing opacities. 

Our next task is to derive the equation for calculating the 
temperature at the end of the timestep at f = t„ [assuming that 
equation ( |35] | has already been solved using some procedure] . 
If we apply some of the approximations in deriving equa- 
tion ( |32] | to the material energy equation dU, we obtain the 
following equation: 



Um,n+\ = Um^n + l-KpAt,, 



kealdfide 



1 JO 



-cf„ kp At„ pU*„ + pSAtn , 



(38) 



where / is the value of / averaged over time interval t„ <t < 
f„+i . The integral on the right-hand-side of equation (l38l l rep- 
resents the total amount of energy absorbed within time in- 
terval t„ <t < f„+i, the 3rd term represent the total energy of 
emitted particles, and the 4th term accounts for the energy ex- 
changed due to physical (not effective) scattering within the 
same interval. All of these quantities can directly be calcu- 
lated by just summing the energies of the emitted and ab- 
sorbed Monte Carlo particles and the amount of energy ex- 
changed in each physical scattering event during a timestep. 

A subtle issue arises here: In the IMC method, effective 
scatterings are introduced in order to model a 1 -/„ fraction of 
the total absorptions and subsequent re-emissions of particles. 
Since the energies of absorbed particles do not necessarily co- 
incide with the energies of the re-emitted radiation particles, 
the effective scatterings should generally be inelastic (as evi- 
dent also from the form of the transport equation l35t . mean- 
ing that radiation particles can exchange energy with material 
due to effective scatterings. However, in equation (|38] | for the 
time update of the internal energy, there is no term that takes 
into account the energy exchange due to effective scatterings. 
Hence, one obvious question to ask is whether it is possi- 
ble to have a consistent Monte Carlo interpretation of equa- 
tions (|35]|-(|38]| if equation ( |38] | does not contain terms that 
account for the energy exchange due to effective scatterings? 
The answer is "yes" if the weights of MCPs are treated in a 
special way during effective scatterings. An approach used in 
IMC photon transport is to assume that the total energy of an 
MCP does not change during an effective scattering, while the 
energy of individual photons within that MCP is allowed to 
change during an effective scattering. Obviously, in this case 
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one has to change the number of photons within that MCP 
in order to conserve the total energy of that MCP during the 
effective scattering. Using this treatment, one can execute a 
consistent Monte Carlo interpretation of equations (l35t -(l38Tl. 
To the best of our knowledge, this feature of the IMC method 
has not been pointed out in the literature previously. 

Once Um,„+\ is obtained using equation (l38T l, the tempera- 
ture r„+i can be calculated by solving iteratively the following 
equation: 



U, 



m,«+l 



C,{T')dT' , 



(39) 



where Cy is the specific heat capacity. Finally, we point 
out that none of the approximations made in deriving equa- 
tion ( |38] | v iolates energy c onservation (see, e.g., FC71 or Sec- 
tion 3.2 of lWollabeiil2008h . 

5.1. Summary of the Monte Carlo procedure. 

The Monte Carlo procedure for solving equation d35i can 
be summarized briefly as follows. Let us assume that / and 
T are known at time t = f„, and we wish to determine them 
at f = f„+i = ti, + At„. The temperature T (as well as other 
relevant quantities such as the opacity k„, etc.) is repre- 
sented in each cell of the spatial computational domain us- 
ing cell-centered and time-centered values. As mentioned 
above, in many practical appUcations the time-centered val- 
ues {Ra,Kp,Ks,>t:s,b,/3} of {Ka,Kp,Ks,:Hs,b,P} in each spa- 
tial zone are given by th e data available at t = t„ [although 
alternatives are possible (IWollaber & Larsenll2009h 1. Using 
these time-centered values, one calculates the sources for each 
spatial zone, generates new particles from the sources, and 
advances both newly-created and census MCPs according to 
the transport equation (l35t by a standard Monte Carlo proce- 
dure (as described in Section |4). In the process of advanc- 
ing MCPs, one keeps track of the total energy of emitted, ab- 
sorbed, and scattered MCPs. At the end of the timestep, {/,„ is 
advanced in time according to equation (l38T l, while the tem- 
perature T is updated using equation i 



6. EXTENSION OF THE FLECK AND CUMMINGS SCHEME TO 
NEUTRINO TRANSPORT 

The FC71 scheme is not directly applicable to neutrino 
transport because, in the latter case, emission and absorption 
of radiation particles not only change T, but can also alter 
the value of Yf. An additional difficulty arises when one has 
to evolve different neutrino types together. In CCSN sim- 
ulations, one usually has to solve three different transport 
equations for three different species of neutrinos: electron 
neutrinos (z^e)> electron anti-neutrinos (P^), and heavy lepton 
neutrinos and antineutrinos, where the latter two are usually 
lumped together into one group (i^J. An additional compli- 
cation arises due to pair processes between neutrinos and an- 
tineutrinos. Such processes couple the two neutrino species 
and add extra non-linearity to the coupling to matter. Here, 
we extend the FC71 equations to the more general case for 
which there are additional degrees of freedom in Y^ and mul- 
tiple neutrino types. For simplicity of illustration, we limit 
ourselves to the case of ID spherically symmetric matter that 
can emit, absorb, and scatter radiation. For such a system, 
the transport equation for neutrinos of type / is again given by 
equation (|3]l, which has to be solved together with the equa- 
tions for the change of the internal energy Um and electron 



fraction Yg: 

P—JT =27ry^ / / Kaiili-Bi)dfide 

+ Y,Si, (40) 

pNA^=2TTY,Si / -^{Ii-Bi)dnde, (41) 

where subscript / is used to denote quantities representing 
neutrinos of type /, and the sum in equations (|40l)-(l4T]) runs 
over all neutrino species. Variable e is again the neutrino en- 
ergy, Na is the Avogadro's number, and i, is a constant equal 
to +1 ,-1 ,0 for i^e, De, and z^^, respectively. Function 5, is the 
function S defined in formula (|5]l for neutrino of type /. 

In order to handle multiple types of neutrinos, we adopt 
an operator-split approach: we evolve different types sepa- 
rately and independently within a timestep. Moreover, we 
linearize the pair processes by assuming that the distribu- 
tion function of the pair counterpart not being followed ex- 
plicitly is given by its local equilibrium value (similarly 
to Janka & Hillebrandt 198 9). This is a good approxima- 
tion, since the pair processes are significant only in high- 
temperature inner regions, where neutrinos are close to ther- 
mal equilibrium JSumivoshi & Yamadal 120 12h . Therefore, 
and hereafter, we focus on solving the transport equation for 
a single neutrino species only. In this case, we will not need 
to sum over neutrino types in equations (I40t-(l4ll: 

dU /■' f°° 

p—^ = 27r / Ka{I-B)dp.ds + S, (42) 

dt J_i Jo 

pNA—^ = 2TTSi / —(I-B)dpde. (43) 

dt J-x Jo £ 

In these equations, we have omitted the subscript /, except in 

Si. 

We now derive some thermodynamic relations which will 
be used later. We first represent the time derivative of the 
specific internal energy U„, using its partial derivatives with 
respect to T and Y^: 



dUm 

dt 



fdU„ 



V dT 



pje 



dT 

'di' 



dU„ 
dYe 



P.T 



dY, 
It 



Using this equation, we can obtain 



dT 

dt 



1 



dU,„ 
dt 



fdU„, 



\dYg 



dYg 

3 y Ul 



^. (44) 



(45) 



where Cy = {dUm/dT^ is the specific heat capacity. Simi- 
larly, we split the time-derivative of U/. 

dt [ dT J „ ^,. dt [ dY, j .,j dt ' 



Using 



and ( |46] |, we obtain the following expression 
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Nher^ 


pCv 


(du,.\ 

KdTj 


p.Ye 


^ pNa 


' (dUr\ 1 


(dU,n 


\ fdU,. 
Jp.AdT 
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p.Ye 



(48) 



(49) 



Note that in numerical simulations, the quantities Cv and 
[dU„,/dYe) can be calculated using the EOSPl The func- 
tion Ur and its partial derivatives with respect to T and Y^, can 
be calculated (semi) analytically using the expression for the 
Fermi-Dirac function (Appendix A). 
We now define two new variables: 

Xa=—, (50) 



and 



Xp-- 



/q°° XaBde 
J-Bds ■ 



(51) 



Following FC71, we again approximate the functions 

{Ka,Kj„K,,>C,,b,Xa,Xp} with {K„,Kp,K„>C,,b,Xa,Xp}^ the 

time-centered values of the former within the time interval 
tn^t < t„+i, and rewrite equations (I42H43) using this approx- 
imation: 

"1 /"OO 

Kaldpde - CKpUr + S , (52) 



dU,„ 

p — : — = /TT 



dt 



1 Jo 



dY, 
pNA-^ = 27rsi 
dt 



1 JO 



xJd^ide-csiXpUr. (53) 



Using (|52] i and ( |53] l and the time-centered values {/3, C,} of 
{/3, C), we rewrite equation (l47T i as 



dUr 
dt 



+c 



27r 



lltSi 



1 JO 

1 



RJd pde — cKpUr + S 



Xaldfide-csiXpUr 



1 ^0 



(54) 



After rearranging terms in this equation, we obtain: 



dUr 



:27r 



{fiKa + (siXa)Idnde 



dt j_, JO 

-0Rp + CsiXp)cUr + ^S. 

For brevity, we now introduce the following notation: 

J=Piia+CSiXa, 
% = PKp+CSiXp, 

and rewrite equation (l55t using this notation: 



dUr 



1 /"OO 



(55) 

(56) 
(57) 

(58) 



■■2tt / ^Id ^de-c^pUr + fiS . 
dt J_i Jq 

" Note the analo gy b etween this quantity and its "photonic" counterpart 
/3 given by equation j28). 

'* In neutrino transport simulations in core-collapse supemovae, the EOS 
in nuclear statistical equilibrium (NSE) is given as a function of three inde- 
pendent quantities (p, T, Ye), usually in tabulated form. 



Next, we apply the time-averaging operator 
tion (|58]l to get 



to equa- 



*~^ r,n+l ^ y.n 



= 27r 



"yldpde-c^pUr + pS. (59) 



Our next task is to eliminate C/,„+i from this equation. In order 
to do this, we make one more approximation. 



Ur = aUr.„+i + ( 1 - a)U*„ , 
which can also be recast as 

U,„^i = U*„ + (U,--U*„)/a, 



(60) 



(61) 



where U*„ = Ur,n + $^t„S and a is the "neutrino" analogue of 
the user-defined parameter a E [0.5, 1] of the Fleck & Cum- 
mings scheme for photons discussed in Section |5] This pa- 
rameter controls the degree of "implicitness" of the method, 
with a = 1 being the most implicit (since Ur = £/r,n+i in this 
case). Substituting C/,.,„+i given by this formula into equa- 
tion ( |59] l, and solving the resulting equation for f7, , we obtain 



where 



Ur = f„U:„+27T 



fn = 



Clp 



1 

1 Jo 



'y Id fide , 



1 + acAtn^p 



(62) 



(63) 



is the neutrino analogue of the "photonic" Fleck factor given 
by formula (|34l i. We now make thefinal approximation of 
FC7 1: we replace the time-averaged Or and /in equation ( |62] | 
with their "instantaneous" counterparts. Or = Ur{t) and / = I(t), 
to obtain: 



t/,(f) = /„t/*„ + 27r 



Clp 



1 nOO 



jl{t)dij,de, (64) 



1 JO 



and rewrite equation (|3]l using the approximation 

{Ka,Kp,K,,>C,,b,Xa,Xp} - {Ka,iip,K,,k,,b,Xa,Xp} dis- 

cussed above. We then have: 

1 dl(p,e) dl(p,e) l-p^ dl{p.,e) 
c dt dr r dfi 

= CKabUr-(Ka+ K!:)I{fl,s) 



+2-K 



ks(e' ,11' — >■ e, p)I(p' ,e')dij.'de' 



(65) 



Substituting equation ( l64l i into equation (l65T l. we obtain the 
transport equation in a new form: 

idi di i-ii^di ^_ ~^, 
■ + /^^ + — :^^r.=f"''^icbUr„ 



dt 



dr r dfj. 

.1 



-(Kfl+Kj)/+27r 



(l-f„)hab 

Ip 



^Idp,de 



I JO 



+27r 



>Cs(e' , p.' -^ e , p)I(p' ,s')dpde' . (66) 



For reasons that will become apparent later, we rewrite this 



equation in a slightly different, but equivalent, form 
\dl dl l-fi-dl 



+ ^ +. 
c at or r 
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We solve this equation for U„i.n+\ and obtain: 
.. Af, 



11 



d^ 



— l^eaCbU i-j^ 



+27r 



kab 



Kes.e^dflde + ln 

1 Jo Xp 



kab 



Xesjidfide 



1 Jo 



/ >fs{£',fJ-' -^ e, fJ-)I{iJ.',e')dfi'de' , 



(67) 

(68) 
(69) 
(70) 
(71) 



where we have introduced a set of new variables: 

0kp 

Kes,e = (1 ~ Jn)~Z ^^a : 

Ip 

a — f \2^jXp 
In) : l^a 7 

Ip 

0kp 

Xes.e — \^ ~ Jn)~Z Xa •> 

Ip 

Xes,/ = (1 -/«)—: Xa, 

Ip 

and Kea is defined as in formula ( |36] |. Equation (|67] |. together 
with boundary and initial (t = t„) conditions for /, determine / 
during the time interval t„<t < t„+i . 

6.1. Update of T and Y^ 

Having derived the transport equation in a new form (|67i . 
our next task is to derive equations for the update of T 
and Ye at the end of timestep t = f„+i, assuming that some 
(Monte Carlo) procedure has been used to solve equation (l6Tt . 
We start by performing the time-averaging integral (l3Ji over 
equation (|52] |: 

p '"'" — ■ = 27r / / kaljdiide-ckpUr + S . (72) 

To conserve energy, we must approximate this equation pre- 
cisely the same way we did in deriving equation ( |67] i. We, 
therefore, substitute equation (l64l into equation (|72] | to ob- 
tain: 

^ 1 />oo 



P ^ — I = 27r 



Af« 



ckp f„U„, + 2TT 



Clp 



1 J(l 

1 />oo 



7/£//^£/£ +5. (73) 



1 JO 



After rearranging some terms on the RHS of this equation 
and using variables defined in formulae (|68)-(|7ll, we obtain 
the following expression 



P — ^ = 27r 

Af« 



1 /'OO 



keJdfids- 



1 JO 
1 poo 



cfnkpUr,n + 2'!r 



kesjldpde 



1 JO 

1 /'OO 



-27r^ / / Xesjidfxde + S. 
Xp J-\ Jo 



(74) 



U,„ 



/ k^Jdjide — 
. , 1 Jo 

/>1 />oo 

cf„kpUr.y, + 271 / / kgsjTdiide 



1 JO 



-27r-^ 



1 /'OO 



Xesjidfide + S 



1 Jo 



(75) 



Using similar arguments, we obtain a similar expression for 
Yy. 



Ye.n+l = Ye.n + — f i llTS, 
RNa 



1 pOO 

/ Xealdfide- 
1 Jo 

1 />oo 

",/„Xpf/r,« + 27ri/ / / XesJdfide 



-2-KSi^ 



1 JO 

1 /"OO 



k^s. eld fide 



1 JO 



(76) 



where Xeo = /nXa- These last two equations determine how 
the values of C/„, and Y^ change after each timestep. 

6.2. Energy and Lepton Number Conservation 

The transport equation (l67t and equations (l75Tl-(l76]l for the 
time evolution of the internal energy U,,, and electron fraction 
Yc conserve the total energy and lepton number in the system. 
This can be demonstrated in the following way. If we apply 
the operator 

^JJ dpde {11) 

to the transport equation (l67l i and add the resulting equation 
to equation ( |75] |. we obtain the following relation 



1 
where 



and 



[Io,n+\ -h,n) + P {Um.n+\ " Umjij 



dh_ 
dx 



Io = 



h = 



Idp,de , 



pidpde 



(78) 



(79) 



(80) 



Clearly, equation (ItHT i is a discretization in time of the follow- 
ing law: 



dt \ c '"I dx 



(81) 



The two terms inside the brackets are the total energy in radi- 
ation and matter, while the term on the RHS is the radiation 
energy flux, meaning that this relation represents the energy 
conservation law. 

Lepton number conservation is also demonstrated in a sim- 
ilar way. If we apply the operator 



1 .OO J 

/ —dfide 
I Jo £ 



(82) 
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to the transport equation (|67] | and add the resuhing equation 
to equation ( |76] |, we obtain a variant of equation (iTSl l for lep- 
ton number (instead of energy), which is a finite-difference 
representation of the conservation law for the lepton number 
Hence, based on this we conclude that none of the approx- 
imations made in deriving the system of equations i&U and 
(l75] )-(l76ll violate energy and lepton number conservation and 
that these two conservation laws are satisfied rigorously. 

6.3. Monte Carlo Interpretation 

We now give a Monte Carlo interpretation for the transport 
equation (|67] | and equations (l75Tl-(l76]l for time evolution of 
the internal energy {/„, and electron fraction Y^, respectively. 

We start with the transport equation (l67t . As in the case of 
photon transport discussed in Section|5] we interpret the terms 
keacbU*^^ and keal on the RHS of equation ( l67b as source 
and sink terms due to emission and effective absorption of 
MCPs. Moreover, terms hes.el and 2tt^ J_^ J^ Kg, Jdiide 
look like terms for a sink and source for scattering. Follow- 
ing FC71, we interpret this scattering as effective scattering. 
Analogously, we assume that the total energies of MCPs are 
conserved in such effective scatterings, while the number of 
leptons in MCPs are allowed to change in order to conserve 
the total energy of the MCP In other words, in such scatter- 
ings, the MCP does not exchange energy with matter, but can 
exchange lepton number 

In addition to these terms, equation ( |67] | contains terms 

Kfs,// and 27r4^ J_^ J^ Xes.i^d iJ.de. These terms again look 

similar to the sink and source terms for effective scatterings, 
but with a subtle difference: In these scatterings, the weight 
of MCPs should be treated differently. Instead of keeping the 
total energy of an MCP fixed, here we fix the total number of 
leptons in the MCPs. This is done in order to be consistent 
with equation JTSl ), as will become apparent in the following. 
Therefore, in these scatterings, the MCPs exchange energy 
with the matter, but not lepton number 

In other words, in order to make a Monte Carlo interpreta- 
tion of equations ( |67l ) and jTSt-jTSl), one has to introduce two 
types of effective scattering. This feature makes this scheme 
slightly different from its counterpart for photons, where one 
introduces just one type of effective scattering. We refer to 
the scattering in which the total energy of MCPs is conserved 
as energy-weight conserving effective scattering, while the 
other type of scattering that conserves lepton number is called 
number-weight conserving effective scattering. 

Let us now consider equations (|75] l and ( |76] l for the up- 
date of the internal energy U,„ and electron fraction Yg, re- 
spectively. Clearly, the 1st and 2nd terms inside the brackets 
on the RHS of equation ( l75b are responsible for the change of 
the internal energy due to absorption and emission of neutri- 
nos within time interval t„ < t < f„+i. Similarly, the 1st and 
2nd terms inside the brackets on the RHS of equation (TEi 
account for the change of Y^ due to absorption and emission 
within the same time interval. Furthermore, the 3rd and 4th 
terms inside the brackets on the RHS of equation ( fTSl l are the 
source and sink terms due to number-weight conserving ef- 
fective scatterings. Analogously, the 3rd and 4th terms inside 
the brackets on the RHS of equation ( |76] ) are the source and 
sink terms due to energy-weight conserving effective scatter- 
ing within t„ < t < r„+i. Finally, the last term on the RHS 
of equation ( fTSl l is responsible for energy exchange due to 
physical scattering, again within t„ < t < f„+i . All of these 



quantities can directly be calculated by summing the energies 
(lepton numbers) of emitted and absorbed Monte Carlo parti- 
cles, and summing the energy (lepton number) exchanged in 
(only effective) scatterings during the timestep. Having calcu- 
lated U,„,„+i and Ye^„+i, r„+i can be obtained via the EOS table 
using the new values of U„,,i,+\ and Yg^„+i. 

6.4. Summary of the Monte Carlo procedure 

The Monte Carlo procedure for solving equation (l67i can 
be summarized as follows. We assume that /, T and Y^ 
are known at time t = t„, and we wish to determine them 
at f = f„+i. The temperature T and electron fraction Y^ (as 
well as other relevant quantities, such as the opacity, k,,, 
etc.) are represented on the spatial computational domain 
using their cell-centered values in each of the spatial zones. 
We assume that the time-centered values {k„, Kp,b,^, . . . } of 
{k,,, Kp,b,'^, ...} for each spatial zone are given in terms of 
the data at the beginning of timestep at f = f„. Using these 
time-centered values, we calculate the sources appropriate to 
each spatial zone, generate new particles from the sources, 
and advance both newly-created and census MCPs according 
to the transport equation (|67t by a Monte Carlo procedure 
similar to the one described in Section |4] In the process of 
advancing MCPs, one keeps track of the total energy and lep- 
ton number of emitted, absorbed, and scattered MCPs. At the 
end of the timestep, we calculate the updated Y,, with equa- 
tion ( |76] |. while {/,„ is updated according to equation (fTsT l. We 
obtain the new value of T using the EOS table with the new 
values of Fg and f/,„. 

7. DISCRETE DIFFUSION SCHEME FOR MULTI-GROUP MONTE 
CARLO NEUTRINO TRANSPORT 

In the IMC method, when the absorption opacity is high, 
the Fleck factor /„ becomes small (/„ ~ 0), and thus k„ — i^a 
and Kga — [cf. equations (|36]|-(|37])1, signifying that most of 
the absorption (and subsequent re-emission) is replaced with 
effective scatterings. In this regime, MCPs undergo Brown- 
ian motion with small (~ I/kcj) mean-free-path most of the 
time. The computational cost of each simulated MCP path 
between collisions is about equally expensive. Thus, simula- 
tions with a large scattering cross section (both effective and 
physical) can be very time consuming due to the large num- 
ber of MCPs paths between scattering events that one has to 
simulate. On the other hand, when the mean-free-path for this 
effective scattering is small, then the solution of the transport 
equation is well approximated by the solution of a diffusion 
equation. Several schemes that aim to make the IMC method 
more computationally efficient at high optical depths by using 
the di ffusion approximation have been sugg ested in the liter- 
ature ( iFleck & Canfieldlll984l: lGentilell2001l) . One of the sim- 
plest and most efficient such met hods is the dis crete-diffusion 
Monte Carlo (DDMC) sc heme of lDensmoreet al. (2007). 

iDensmore et al.l (l2007h developed the DDMC scheme for 
gray radiation transport without physical scattering in ID pla- 
nar geometry for non-moving matter In this section, we ex- 
tend this scheme to the energy-dependent case with physical 
scattering (the extension to the velocity-dependent case pre- 
sented in Section |8]l. We again assume a ID spherical static 
matter distribution and, for simplicity of illustration, focus on 
photon transport (instead of neutrino transport) because the 
ideas behind extension to the energy-dependent case do not 
depend on any aspects that are specific to photons or neutri- 
nos. In the following, we first derive the discretized diffusion 
equations in the multi-energy case with physical scattering. 
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and then give a Monte Carlo interpretation of the relevant dif- 
fusion equations. 

We start by introducing the zeroth and first radiation mo- 
ments J and H (.Mihalas & Mihalas 1984.) : 



/=- / W/i, 



H=- 



and apply the operator 



1 



I^jid^ 



dn 



to the IMC photon transport equation 
IdJ Id 



(83) 
(84) 

(85) 



to obtain 



■ + ^ — [r H) = fnKabcU*„ - {ka + Ks)J 

c at r-^ or ^ ' 



+4^(1 -/«)^ / lJ(e')de' 
Ip Jo 



+7r 



+1 /.+1 <.oo 



i<i(e',/i' -^ e, ^)I(^',e')d^d^'de' .(86) 



I JO 



Observing that 



2tt 



+1 



— .-,0/^/ 



ksie ,;U — > £, ^J,)d^ = k^ie — > e) 



(87) 



does not depend on the "angle" /i, we rewrite equation 
using this result: 

I gj I g 

c at r^ or ' 



+47r(l-/„)^ / lJie')de' 

+ / x°(e' ^ e)J{e')de' . 
Jo 



(88) 



We assume that our entire spatial domain < r </? is di- 
vided into connected, non-overlapping spatial cells, for which 
r G ['";-i/2i 0+1/2]' where j = {!,... ,A^;.}. Furthermore, we 
designate a subregion < r < R^^ of this domain, which is 
covered by y = {1, . . . ,m} cells, for DDMC. The cells with 
j < m will be called interior DDMC cells, while cell j = m 
will be called the interface DDMC cell. The discretized dif- 
fusion equations for interior cells are slightly different from 
those for interface cells. Therefore, we derive them in two 
separate steps in the following. 



7.1. Interior DDMC Cells 

For interior DDMC cells, we approximate equation 
in each spatial cell j by using cell-centered values 
{Ra,j,Ksj,f„j,U;^j,bj,-/j,-/pj,>i:lj} of the quantities 



{K.a,K,,f„,U*„,b,J,Jp,3^}: 

\_dJ_ ^_9 
c dt r^ dr 



I dJ I d 

+ — T^ (r^H) = fnRa.jb cU*„ j- (kaj + ksj)J 



+4TT(l-f„j)^''^' 



Ip.j Jo 



^jj{e')de' 



+ 
Next, we apply the operator 

dV: 



/ k^ j{e' -^ e)J{e')de' . (89) 

JO 



j+1/2 



r dr. 



AVj 7,^_,^, E.Arjrj J,. 

where Ej = 1 + Ar?/(12rJ), to equation ( [89] ) to obtain 
1 dJj 1 / 2 2 \ 

c^^ s;a7;;j v>'/2^i+i/2-o-i/2^i-i/2J 

— Jn l^a.jb CU i-jjj — (KaJ + Ksj)Jj 
K b- f°° 

+4^(l-/«,>)4^ / ljJjie')de' 



(90) 



k' '-' 



where 



and 



Jj = ^ 



^An>j.r,_,, 



Hj±i/2 = H{rj±i/2)- 



Ip.j Jo 
j{e' -> e)Jjie')de' , 

r^Jdr 



(91) 



(92) 



(93) 



We further transform equation (l9Ti for the interior DDMC 
cells. Using Fick's law jPomraningll 19731) . 



H(r): 



1 dJ 



(94) 



3kt dr ' 

where kj is the transport opacitjQ we evaluate the cell-edge 
Hj+i/2 of H within time interval t^ <t < f„+i as: 



H 



J+l/2 '■ 



-5^a;('--OW2). 



(95) 



while Hj_ii2 is calculated similarly at point r = ?"y_i/2- Note 
that here we have used the time-centered value kt of the trans- 
port opacity kt- By employing a finite-difference approxima- 
tion for equation ( |95l l, we can express ^;+i/2 in cell j as 



H 



j+^12 '■ 



^^Tj-H/2^0 



or in cell j+l as 



H 



j+1/2 - ■ 



3k 



TJ-+1/2 



Ao+, 



{jj+l/2-Jj) 



{jj+l-Jj+\/2) 



(96) 



(97) 



where Jj+1/2 = J{f = 0+1/2) ^nd k^ .^j .j is the transport opacity 
at the inner boundary of cell 7 + 1 , while k 



T,;+l/2 



is that at the 



' The transport opacity kj is defined as kt = Kc + Ks-Itt J^ fj,Ks{fi)d^. 
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outer boundary of cell j. Equating the RHSs of equations 

and Wi\ . and solving the resulting equation for Jj+in, we ob 
tain 
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RHS of this equation as 



J 






(98) 



J+l/2 '■ 



Then, if we use equation (|98l l to evaluate either equation i 

or equation ( |97] i. we find an approximate expression for 



4^(1 -/;,j) 

■■47r(l- f„,j) 



Kaj(Sk)bj(ek) 



^^j(ei)Jj(ei)Aei 



IpJ 



^^j(ei)Jj(ei)Aei 



't* 



H 



J+l/2- 



Hj+i/2 ■■ 



Jj+l-Jj 



+(1 -fn,j)Haj(ek)Jj(ek)-^ — /.-.., , (104) 



(99) 



and 



Substituting the RHS of ( |99l ) and a similar expression for 
Hj-i/2 into equation ( |9T] l. we obtain an equation for Jj in cell 

l^y,. = -[.,, + .,,+R„, + ;i,j7,.+/„,S„,Vf/:,, Y.~<M'^'^'^-'Aei)^e, + >cl,{ek^ek)Jj{ek)Aek. 



/o 7;(£)5j(e¥e^ 
^>;r°/£,^ei)7/£,)Ae,= 



dt 



l¥k 



-M^rj+yr]+i ^j.yArj.\r]_^ 



(105) 



^Anr] 



-KL,j+lJj+l-i — 



^An>-] 



KR,j-lJj- 



+47r(l-/„j)4 



'^"'^^^ I jjJj{e')de'+ I k'^j{e'-^e)Jj(e')de' 
Ipj Jo Jo 



' Using equations ( I104b -( I105I ) in equation ( 1103b . we obtain 

Jj,k = -[KLj,k + KRj^k + fn.jK.aJ,k + {^ -fn.j)^a.jM 



(100) 



In the last equation, we have introduced two new quantities: 

1 



dt 



S,-„Ar--2 



+S'sJ,k] Jj.k+fnjiia.j.kbj,kCU*„ J 






"I =—. =; I^L,j+\,kJ i+\M T ^ 2 ^R.j-lMJj-l.k 



and 



'*«,, 



35;Aor2 li^._^^^Arj + R-._^^^Arj.i ' 



Hi/2 ^ 

3^jArjrj R-j^^^^Arj + R^j^^^^Arj+i ' 



(101) 



(102) 



-j^'j'j 



^aj,kt>j^k ' 



+47r(l-/,,,)^^2^^^^^^.,y.,A„ 



7p, 



/# 



which are called left-leakage (kij) and right-leakage (krj) 
opacities. The reason they are called "leakage" opacities will 
become apparent below, when we provide a Monte Carlo in- 
terpretation of equation (llOOt 



(106) 
where the subscript k denotes quantities pertaining to energy 



As a next step, we discretize equation ( fTOOl i into energy group k. In equation ([T06l), for brevity we also have intro- 



groups: 
1 d 

- JfJA^k) = - [l^L,Mk) + KRjiSk) + Ka.j{£k) + R,j(ek)] J jiSk) 

+fn.jiio.jiek)bjiek)cU*„ . + — ^—- — ^ — KLj+i(et)/y+i(£/i) 

^jArjrj 

Ej-iArj-irj_i , ,, , , 

2 I^R.j-l(£k)Jj-l(£k) 



duced two new variables: 



and 



CTaJM = 



(Ts.jM = 



^ jj.kBj,kAek 
J^-fj{e)Bjie)de 



>t^ :{ek^ek)Aek 
1 



K.j.k (107) 



^•'<,i,k ■ 



(108) 



+47r( !-/„,,) 



K.aj(ek)bj(ek) 
lp,i 

+ H^°/^' ^ £k)JMi)Ae, 



It is easy to see that da.j.k < ^^a.j.k and cFs.j,k < f^sjM- As we 

will discuss below, this property has important implications 

y^7 (e;)7 (e/)A£; for 'he computational efficiency of our multi-group DDMC 

, scheme. 

We observe that equation (1106b can be viewed as an equa- 
tion for the time evolution of y^jt incell j and energy group k. 
Namely, according to equation ( |106b . function Jjk decreases 
at a rate 
(103) ^ ^ 

...^ . ,JA-u-JU [KL.j.k + KR,j,k+fn.jKa.j.k + (l-fn,j)<Ja.j.k + (T,.j.k\Jj,kC, (109) 

where e; is the value of energy m group / and Ae/ is the width 

of that energy group. We now express the summations on the due to the 1st term on the RHS of that equation and increases 



Monte Carlo Neutrino Transport 



15 



at a rate 



cally using a formula similar to equation 



Jn,jHa,j,kt>j,kCUi.„ j+ ^— 2 '^LJ+\,kJj+lM + 



S;_iAo_ir2. 



EjArjPj 



■I^RJ-\,kJj-\,k 



+4^(1 -/«j) 



KaJ,kbj,k \-^ ~ I \ 



l^ 



lik 



.-,0 



■,./£/ -^ £k)JjjAei 



c,(110) 



due to the rest of the terms on the RHS of equation ( llOOI i. 
Now, recalling that function Jj^k represents the number 
of MC particles in cell j in energy group k, we make 
the following Monte Carlo interpretation of this equation: 

The terms f„jka,j.kbj.,kC^U*„j, 

Sj_iAr,_ii-J_| 



T^oT 



-c^l,j+\mJj+\ 



and 



-cKRj-i^kJj-i describe the rate of increase of the 

number of MCPs in cell j and energy group k due to emis- 
sion and leakage from the right and left neighboring cells, 
respectively. Moreover, terms fnjUaj.kJj.kC, nijj^JjkC and 
KR,j,kJi,kC represent the rate of decrease of Jjk due to absorp- 
tion of MCPs, and leakage of MCPs to left and right neigh- 
boring cells, respectively. In addition to these terms, we have 
terms that are responsible for physical and effective scatter- 
ing. Terms (7sj,kJj.k and (l- fi,.j)o'aj,kJj.k describe the de- 
crease of y^jt due to physical and effective scattering of MCPs 
from energy group k to any other energy group / (l ^ k), re- 
spectively. Accordingly, the terms 

4^( 1 _ /„ .) I^thl. J2 lj,kJjJ^Ae, (111) 

'^P'j lik 

and 

^^°/e,^ £*)/;,, Ae/c (112) 

/^ 

represent the increase rate of J j_k due to scattering of an MCP 
from energy group l^kio energy group k. 

In this picture, DDMC particles have no propagation angle 
^ or position r within a cell, but they always know their cur- 
rent cell and time. A DDMC particle can either remain in its 
cell within a timestep without a collision, or undergo a "colli- 
sion." Here the term collision refers to an absorption, a leak- 
age to the left or right neighboring cell, or from one energy 
group to another The quantity 



l^m.k + I^R,j,k + fnjfia,j.k + ( 1 " fn,j)^a,j,k + &sj,k 



(113) 



can be regarded as the total collision opacity. Using this in- 
terpretation, we can perform DDMC transport based on the 
calculation of distances, similar to the Monte Carlo procedure 
described in Section |4] However, in DDMC, we do not cal- 
culate distances to the boundaries. Instead, we calculate the 
distance to collision, dc, and distance traveled until the end of 
the timestep , d,. Since the distance to collision dc is based on 
the opacity ( 1113b . this distance can be calculated probabiUsti- 



dc = - 



In^ 



'^L,i,k + K,R,j,,k + fn,i'i^a,i„k + ( 1 " fn.j)^a,j,k + ^sj,k 



(114) 

where ^ is a random number uniformly distributed in (0,1], 
and the distance traveled to the end of the timestep is calcu- 
lated using relation ( |19ll'^'1 . 

If the time to collision is less than the time remaining in 
the timestep, the DDMC particle undergoes a collision, and 
the time to collision is decremented from the time remaining 
in the timestep. Again, as w e see from the second term on 
the left side of equation ( 1 1001 ). a "collision" can be an absorp- 
tion, a left-leakage, a right-leakage, or a leakage from one 
energy group to another The collision type is sampled from 
the probability of the collision type that is calculated using the 
relative magnitudes of the different "opacities." For example, 
the probability of left-leakage can be calculated from 



PL- 



KL.j,k 



HLJ,k + KRJ.k+fn,jKa,j,k + ( 1 " fn.j)^a,j.k + (ysj,k 



(115) 



If the collision is an absorption, the MCP history is termi- 
nated, as in standard Monte Carlo. If the DDMC particle un- 
dergoes a leakage reaction, it is transferred to the appropriate 
neighboring cell, and the simulation continues. If the colli- 
sion type is the "leakage" from one energy group to another, 
then the new neutrino energy is sampled using the functional 
form of the differential scattering opacity. If the time to col- 
lision is greater than the time remaining in the timestep, the 
DDMC particle reaches the end of the timestep and is stored 
for simulation in the next timestep. 

We point out that the DDMC a ppro ach is based on the diffu- 
sion approximation to equation (|35] |. so it should yield accu- 
rate solutions when used in optically-thick regions. As men- 
tioned above, the DDMC transport process consists of discrete 
steps that reflect transfer of MC particles between spatial cells 
(but not between spatial locations within a cell, as in a pure 
MC method). Due to this property, DDMC can be much more 
computationally efficient than the standard Monte Carlo im- 
plementation of equation dSSl ). 

It is interesting to note that if the physical scattering is 
elasti c, th en in-scattering and out-scattering terms in equa- 
tion dlOOl ) cancel each other Hence, the presence of elastic 
scattering does not lead to the appearance of any new terms in 
equation dlOOb compared to the case when there is no physi- 
cal scattering. Instead, the scattering modifies only the values 
of the leakage opacities kJ and kJ which, in turn, are a 
result of the use of the transport opacity in Pick's law. For 
this reason, the DDMC scheme does not need to perform any 
special explicit numerical operation in order to model elastic 
scattering; the effect of elastic scattering is taken into account 
via modification of the values of the leakage opacities. For 
this reason, the DDMC scheme leads to the biggest savings 
in computational cost in regimes dominated by elastic scatter- 
ing. 

There is also another reason for higher efficiency of DDMC 
compared to MC schemes, which stems from the following. 
In a regime where the absorption opacity is high, the effective 

-" This procedure is sliglitly different in the continuous absorption method. 
In this case, we again calculate two distances. However, the distance to colli- 
sion is calculated using the opacity k,l.j + f^Rj + i^ ~ fii.j)^a,j + Hs.j and, thus, 
a colUsion can be either left-leakage, right-leakage, or physical or effective 
scatteri ng, w hile the condition for absorption is calculated as described in 
SectionKj] 
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scattering opacity a-^s./.k = {^-fn.j)^a.j.k will do minat e the col- 
lision opacity for DDMC given by expression (I113l l. Hence, 
the effective mean-free-path in this regime is ~ I/ctcj.j.a- 
However, recalling that aaj.k < iiaj,k, we see that the effec- 
tive mean-free-path for DDMC should be larger compared to 
that for the MC scheme, which is given by 1/[(1 -fnj)iia.j.k]- 
Depending on the energy group k and its width, 5-„ jjt can be 
significantly smaller than Kajk- Therefore, from this source 
alone, we gain a speed-up by a factor of < iia.j.k/^a.j.k- Simi- 
larly, we achieve speed up from the fact that the inelastic phys- 
ical scattering opacity in the DDMC regime, dsj^/;, is smaller 
than that for the MC scheme, which is ftsj-,*. 

Finally, the biggest speed-up in DDMC comes from the fol- 
lowing assumption at a cost of making one more (but excel- 
lent, as will become clear later) approximation. We split the 
effective scattering opacity d'esj,k into two parts, Uj^kO'es.j.k and 
il-aj.k)d'esj,k, where < Qj^i^ < 1, and restrict the first of 
these two to be elastic effective scattering, while the second 
one is free to be inelastic (as effective scattering would oth- 
erwise be). As we discussed in the above, the presence of 
an extra elastic scattering source does not increase the cost 
of doing DDMC transport. Therefore, by assuming that a 
Oj,k^es.j.k fraction of effective scattering is elastic (instead of 
being inelastic), we achieve computational savings propor- 
tional to aj,k- As we demonstrate in Section W2\ depending 
on the scenario, this can lead to speed-up of calculations by a 
factor of 10^- 10"^ or more. 

Since effective scattering is in general inelastic, an obvious 
question arises: Is it a good approximation to treat a fraction 
cfyA of effective scatterings as elastic in the DDMC region? 
The answer depends on the value the Fleck factor f„j and the 
parameter fl J jt. The only way the inelasticity of effective scat- 
tering affects the transport of MCPs is by enabling thermal- 
ization of MCPs when they move to a new cell with different 
T, Ye, or p. By thermalization we mean a change of the spec- 
trum of MCPs when they move to a new cell to reflect the 
emissivity spectrum of that cell. However, when each MCP 
undergoes inelastic scattering at least once after it moves to a 
new cell, that MCP acquires an energy spectrum that reflects 
the emissivity spectrum of its new cell. In this case, the inelas- 
tic nature of effective scattering should not play any role after 
MCPs move to a new cell and before they leak out to another 
zone. Thus, these scatterings can be treated as elastic between 
the two events. Alternatively, if MCPs that move to a new cell 
get absorbed before they propagate to another zone, the in- 
elastic nature of effective scattering again should not play a 
role. Therefore, for this treatment to be exact, MCPs that leak 
to a new cell should undergo effective scattering at least once, 
or get absorbed before they move to a different cell. This con- 
dition is fulfilled if 



or 



K-L,j,k<^fn.j&a,j,k and 



'^R.j,k-^fnJ<^aJ,k, 



KLj.k < aj,k(l - fnj)a„j,k and 



(116) 



7.2. Interface Cells 

The method which we u se for interfacing DDMC with stan- 
dard MC is the same as in lDensmore et al.( (120071) . 

We start by deriving an equation for the cell-centered value 
of Ji„ in cell m on the right boundary of the DDMC region. 
Using equation (|99] l, we derive an expression for H„j-i/2'- 



H, 



^m «'; 



m-1 



m-1/2 - " 



(118) 



Substituting ( |118t into equation ([91) for cell j = m, we obtain 



£ CA/m __r -1. 

C Ot '- J „ 



' l+m/2 



^m^^m^m 



H, 



m+l/2 



/■ ~ 7 TT* , '—'11-1 ^ ''m-1 ''m-1 T 

'T'Jnl^a.mOCU i.j^ ij^-\ — - ^ KR.m-lJm-l 



~ Ar r~ 



47r(l-/„.J-:^ / 7,„7„,(£)'i£ 

Ipjn Jo 
/•oo 

+ / k'l„(e' -^ s)Us')de' Xim 
Jo 

where we have made use of equations dlOlb and (I102l i. To 
complete this derivation, we must find an approximate expres- 
sion for the flux H at the in terface of the DDMC region. 

Following iDensmore et a l. (2007), we use the asymptotic 
diffus ion-limit boundary condition dHabetler & Matkowskvl 
fT975h : 



1 X 

W{p)Ib{p., t)dfi = J(r,„+i/2) + : 



df 

dr 



r=r„Mll 



'0 "'T.m+1/2 

(120) 

where hiii^t) is the radiation intensity due to Monte Carlo 
particles incident on the DDMC region, A ~ 0.7104 is a con- 
stant, and W(/i) is a transcendental function well approxi- 
mated by 



(121) 



Incident intensity Ii, in equation (1120b is weighted by W{p), 
which takes into account the angular distribution of the MC 
particles coming into the DDMC regi on. 

To express //,„+i/2 using equation ( 11201 ). we approximate 
the derivative on the RHS of equation ( 11201 ) with a finite dif- 
ference: 

/•I 2A 
/ W{p)Ib{p.,t)d^i = J„+ii2 + - T — (4,+i/2--/m), (122) 

where kjm is the cell-averaged value of the transport opacity 
k-Y in cell m, and Jm+\/i is an appr opriately defined cell-edge 
value of y. Solving equation ( 11221 ) for /„,+i/2' we obtain 



J, 



m+l/2 ■ 



KT,,„Ar„, + 2A7o 



W(/i)4(/i)£//i+ — 



2A 



~Jm 



kjj„Ar,„ + 2X 

(123) 
Next, we use equation ( |96] | to represent H„,+i/2'- 



KRj,k < aj,kil -f„j)aaj,k. 



(117) 



In Section 19721 we demonstrate that these conditions are met 
in DDMC regions for appropriately chosen values of ay /.. 



H, 



m+l/2 ■ 



3KT.n,Ar„ 



J, 



m+l/2' 



Jm 



(124) 



where we again have used the cell-averaged value kt,,, of the 
transport opacity Rj. Substituting the RHS of equation ( 11231 ) 
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into formula (1124b . we obtain 



H, 



m+l/2 ■ 



3KT,mAr,„+2A \Jq 



W(fi)Ufi)dfi + J„, . (125) 



Substituting the RHS of the last equation into equation ( |119t . 
we obtain 

i OJjfi r ™ ™ "I J 

c at '- -' 

'Jn.m^a.m'^m ^^ r.n,m ^ a 2 ^R^m—lJm—l 



'm+\/2 



- / PQi)^iIb(^J.)dn 



^m^riiir^j Jo 



+47r(l-/„,„,)^^ / T^/KeVe' 



Ip.m Jo 



+ / >r?,„(e'->£)^(eV£', (126) 
Jo 



where the right-leakage opacity is defined as 

1 



kr. 



2''„,+ l/2 



^m^^m^m T m+l/2 '^^ 



instead of equation ( 11021 ). and P(fi) is defined as 

2 



PW-- 



T,m+l/2 '"~'~ ^ 



, 3 
1 + -^ 



(127) 



(128) 



Finally, we discretize equation (I126l l in energy groups to ob- 
tain 



1 dJ,n 



c dt 



■ \j^L,m,k + f^R,m,k + l^a,m,k + '^j.m.A-J Jm.i 



'Jn^m^a^m^k'^m^k^^rnm ^ 



^m-1 ^^m-l ^m-1 



" Ar r2 
' — ui'—^'m'tn 



'm+l/2 



I^R,m-l,kJm-l,k 



- / Pk{fJ')f^h,k{fJ-)dn 



+47r(l -/„,„) — 2_^ j,„jJ,„jAe, 



lp,m 



't* 



(129) 

where subscript k is again used to denote quantities pertaining 
to energy gr oup k . 

Equation ( 11291 ) is similar to equation ( 1106b . The only dif- 
ferences are the expression for the right-leakage opacity KR.„,jt 
and the presence of the source due to MC particles coming 
into the DDMC region. The flow of energy due to this in- 
coming radiation for a direction /i is given by j^di,. Therefore, 
P(/i) can be interpreted as the probability with which an inci- 
dent MCP w ith direction ^i converts into a DDMC particle. 

Following iDensmore et al.l ( 120071) . we implement the con- 
version of the MCPs into DDMC particles (and vice-versa) in 



two separate ways, depending on whether the DDMC bound- 
ary is at the problem boundary or no t. In the latter case, we 
use the probability given by equation (1128b to determine if the 
incoming MCP is converted into a DDMC particle. If con- 
verted, it starts transporting using DDMC in cell j = m. Oth- 
erwise, the particle returns isotropically to cell j = m+\. The 
DDMC particles that undergo right-leakage reactions from 
cell j = m to cell j = m+\ are also placed isotropically at the 
boundary of the DDMC region (i.e., at the inner boundary of 
the cell j = m+ 1). Note that this angular distribution is correct 
only when the incident intensity is nearly isotropic. Hence, it 
is important to choose the boundary between the DDMC and 
MC regions where the distribution is sufficiently isotropic. 

Second, if the DDMC region is at the outer boundary of the 
system, then the incoming MCPs are regarded as a particle 
source due to boundary conditions. In this c ase, we split in- 
coming MC particles according to equation (1128b : a fraction 
P(/i) of these particles is converted into DDMC particles and 
begins transporting using DDMC in the DDMC region, while 
the remaining fraction 1 -P(/i) is regarded as MC particles 
escaping the system. 

8. VELOCITY-DEPENDENT MONTE CARLO AND DDMC 

Thus far we have discussed radiative transfer in material 
that is not moving. In this section, we extend the schemes 
discussed in the previous section to the case when matter is 
moving with an arbitrary velocity. 

Again, we assume a spherically symmetric distribution 
of matter, split our spatial computational domain into non- 
overlapping zones, and linearize the transport equations 
within a timestep t„ < t < tn+i- Moreover, for simplicity of 
illustration, we consider the case with no physical scattering. 
However, as will become clear in the following, inclusion of 
scattering is conceptually simple and our code is capable of 
handling physical scattering. We also focus on photons be- 
cause the ideas behind the extension to the velocity-dependent 
case is the same for both photons and neutrinos. We assume 
that the radial component V,. ^ of the velocity vector V (as well 
as other information such as temperature, density, etc.) in 
each cell j does not change within a transport timestep, while 
the and (j) components of the velocity are assumed to be 
zero everywhere. Since MC and DDMC methods are based 
on somewhat different techniques, we separately discuss the 
extension of each of these to the velocity-dependent case. 

8.1. Velocity-dependent MC scheme 

Our velocity-dependent MC scheme is b a sed o n the 
mixed-frame formalism of iMihalas & KleinI ( 119821) and 
iHubenv & Burrowsl (l2007h . In this formaUsm, emissivities 
and opacities are defined in a frame comoving with the fluid, 
which are then Lorentz-transformed to the Eulerian lab frame, 
in which transport is performed. 

Before we describe our velocity-dependent MC algorithm, 
we present formulae for the Lorentz transformation between 
the comoving and the lab frames for several quantities that 
will be useful later in the section. The four-momentum of a 
massless radiation particle is given by 



M" = -(l,n), 
c 



(130) 



where e is the photon (or neutrino) energy, and n is 
a unit spatial 3-vector in the particle propagation direc- 
tion. In spherical polar coordinates, M" has the following 
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form ( iMaalas & Maalasll 19841) 



M° 



l,^,(l_^2y/2^^(l_^2y/2^ 



rsinf 



(131) 



where Lp is the azimuthal angle. If a particle has energy e and 
travels in direction {/^t, ip] as measured in the lab frame, it will 
have some other energy Eq and direction {/io, "^o} as measured 
by an observer attached to a fluid element moving with veloc- 
ity vector V relative to the lab frame. (Hereafter, we denote 
all of the quantities measured in the comoving frame with sub- 
script 0.) Because M" is a four- vector, its components in the 
two frames moving with respect to each other with velocity V 
are related by general Lorentz transformations. Therefore, we 
obtain 

£o = 7£('l-— V (132) 



and 



no = 



eo 



-7- 



7 n- V 

7+ 1 c 



(133) 



where 7 = (1 -y^/c^)~2 is the Lorentz factor In the spher- 
ically symmetric case, for the special case of matter motion 
along the radial direction, these relations take the following 
form; 

eo = 7£fl-— V (134) 



l-^V,./c 
(^0 = 1^, (136) 

where V^ is the radial component of the velocity. Next, we 
need a formula for the transformation of opacity n (both for 
scattering and absorption), derived first by iThomasI (119301 ): 



k(M:£)=— '«o(eo) ■ 



(137) 



From the information about the fluid in each cell (such as 
temperature), we calculate the emissivities and opacities of 
the material in that cell as measured by an observer comov- 
ing with the fluid. We then calculate how many particles are 
emitted in the comoving frame in each cell, and sample the 
radial coordinates ro, propagation directions /^o, energies eo, 
and emission times fo of these newly emitted MCPs in the co- 
moving frame using the comoving frame emissivities in same 
way we did in Section |4] Next, we transform quantities £0, 
/^io, A), fo for each MCP to the lab frame, where we then 
transport the M CPs. The particle energy is transformed using 
form ula (I132l i, while the angle [i is transformed using equa- 
tion ( 11331 ). In order to transform the radial coordinate ro, we 
assume that, at the beginning of a timestep f = f„, the radial co- 
ordinates of the inner boundaries of the comoving and the lab 
frame cells coincide with each other Then, the radial location 
ro of an MCP in the comoving frame in cell j is related to the 
lab-frame radial coordinate r via the Lorentz transformation: 



^7j['-o + V,,j(fo-f„)] , 



(138) 



where V^j is again the radial velocity of the fluid in cell j 
measured in the lab frame and 7^ is the Lorentz factor in zone 
j. The MCP emission time is transformed into the lab frame 
using the formula 



f = 7y h-tn + 



Vrjro 



(139) 



Now, having transformed all the necessary information 
about particles into the lab frame, the next step is to transport 
the particles in this frame. The transport algorithm is similar 
to that for static matter described in Section |4] and is again 
based on the calculation of the distances to collision, spatial, 
or time boundaries. However, in this case, the distances to 
collision (absorption or scattering) need to be calculated us- 
ing the lab-frame opacities, which are calculated from their 
comoving frame values using formula (1137) . If the smallest 
of the distances is d,, then the MCP goes into the census for 
the next timestep, as in the static case (Section|4|. If di, is the 
smallest distance, then the MCP moves to the new cell, where 
we transport the particle using the lab-frame opacity of the 
new cell. If the MCP is absorbed, then its energy and momen- 
tum (and lepton number, if we are dealing with neutrinos with 
lepton number) are deposited into its current cell. If an MCP 
undergoes scattering, then we first transform the energy e and 
angle /i of the particle into the comoving frame, calculate the 
new values of eo and /io as a result of scattering (as described 
in SectionlH, record how much energy and momentum is ex- 
changed between the MCP and matter as a result of scattering 
(in order to deposit both in that cell at the end of the current 
timestep), and transform the new energy eo and angle juo back 
to the lab frame to continue the transport of the MCP. 

8.2. The velocity-dependent DDMC scheme 

Since it is most natural to formulate the diffusion equations 
in the Lagrangian frame, we perform velocity-dependent dis- 
crete diffusion Monte Carlo transport in the comoving frame . 
We start with the equation for the energy density ErPl of the 
specific intensity in t he comoving frame accurat e to 0{Vr/c) 
[cf. equation 95.82 of lMihalas & Mihalasl(ll984l) 1: 

Dt oM,- r Dt 



_d_ 
den 



eo 



r Dt 



= ko(4ttB-cEo),{UQ) 

where D/Dt is the Lagrangian time derivative, d/dM,. = 
l/{4TTr^p)d/dr, Fq is the radiation flux {Fq = A-ttHq, where Hq 
is defined as in equation ([84])), Pq is the diagonal component 
of the radiation pressure tensor. 



1 



In'n^dn, 



4ir 



and 



Dlnp ^ 1 d , 2y\ 
Dt r2 dr ^^ '' 



(141) 



(142) 



Due to the DDMC approach, we assume isotropy of the ra- 
diation in the comoving frame, which impli es tha t Po = l/3£'o- 
If we substitute 1 /3£'o for Po in equation ( I14QI ) and drop all 
terms of 0{XpV,-/lc) and higher, where Xp is the mean-free- 
path and I is the problem domain size, we obtain an equation 
for the evolution of £0 that is valid in the non-equilibrium dif- 

-' Note that Eq = 4ttJq/c, where Jq is defined as in equation (83). 
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fusion limit for moving matter (iMihalas & Mihalaslll984l) : 

d 



D_Eo l_ 
Dt p 3 

1 d 



Bo- 



de 



-{£()£()) 



D 1 



+^^(rVo) = Ko(4^B-c£o) 



(143) 



It is easy to show that the last equation can be cast in the 
following form: 



DEq 
Dt 



^D\np eodEoDlnp 1 d ,2r^\ 



■■ko{4ttB-cEo). (144) 



Using equation ( 1142b and recalling that D/Dt = d/dt + V,- ■ V, 
we rewrite the last equation as 

dEj) dEo 9\^ ep dEpDlnp 
dt 'dr °9r 3 9eo Dt 



I d , ^ 



FQ)=Koi4nB-cEo). (145) 



This equation incorporates three different velocity- 
dependent effects: advection, number density compres- 
sion/decompression, and Doppler shift, which are described 
by the 2n d, 3r d, and 4th terms, respectively, on the LHS of 

equation (1145b . 

We now compare equation (1145b with the other equations 
for isn (under the same approximation) used in the literature. 
ICastod (120041) derived an equation for £0 accurate to 0{V /c): 



f^ + V.(V£o) + 



P(y 



d(ePo) 
de 



: VV + V-Fn 



= 47r;o-Kac£'o , (146) 



where Jq = K()B and the colon ":" operator indicates summing 
the product of the tensor on the left with the tensor on the right 
over two indices, viz.. 



R:S 



'■ Z^^ijSij ■ 
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Rewriting equation (1 146b in ID spherical polar coordinates 
and substituting Po = l/3£'o, which is valid in the diffu- 
sion limit , we obtain an equation that is equivalent to equa- 
tion ( 1145b . 

Swestv & Mvral (120091) derived an equation for £0 accurate 
to 0{V/c): 



dEo 
dt 



+ V-(V£o) + V-Fo 



djPo : VV) „ 

'^ 5 = ^ 

ae 
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where S is the collision term. If we substitute Pq = l/S/io into 
this equation and rewrite the resulting equation in ID spher- 
ical polar coordinates, we obtain an equation equivalent to 
equation ( |145b . 

8.2.1. Numerical Implementation 

To implement the DDMC scheme, we first rewrite equa- 
tion ( 1145b using the quantities J = cE/{4tt) and H = F/{4tt) 



that we used in Section |7] 



1 dJo Vr dJo Jo dVr £0 dJo Dlnp 
c dt c dr c dr 3c deo Dt 

+^^{rHo)=Ko(B-Jo). (149) 

r~ dr ^ ' 

It is easy to show that an IMC version of this equation has the 
following form: 



1970 y,dJQ JodVr eodJoDlnp 1 
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We now split this equation into three separate equations: 
d 



I d Jo j__ 
c dt r~ dr 
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and 
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1 = 0, 

c dt 3c deo Dt 



IdJo^V^dJ^^JodV^^Q 
c dt c dr c dr 



(151) 



(152) 



(153) 



and solve these three equatio ns in operator-split manner in 
three separate steps. Equation ( 1151) is responsible for the dif- 
fusion of radiation through matter and is the same as equa- 
tion ( l86b for DDMC transport for non-moving matte r He nce, 
in the first operator-split step, we solve equation (1151b the 
same way we did equation ( [86] ). 

Equation ( 11521 ) is responsible for the Doppler shift of the ra- 
diation energy. In order to solve it, we first rewrite this equa- 
tion using a new variable eo = Ineo: 



dJo dJo _ 
dt deo 



where 



= 



IDlnp 
3 Dt 



(154) 
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Equation ( 1154b is akin to an advection equation. Its solution 
within a timestep (within which g is assumed constant) at time 
t is given by 



Joit) = J'^\eo-git-to)), 



(156) 



where /™ is the value of Jo at initial time fo. Due to this sim- 
ple analytical nature of this solution, there is an easy way of 
incorporating equation (1154b into the full velocity-dependent 
DDMC framework. Whenever we perform a DDMC trans- 
port operation on an MCP (i.e., move an MCP from one spa- 
tial cell or energy group to another) within a time interval At, 
we shift the energy of the neutrinos (or photons) in that MCP 
by an amount corresponding to Ae^ = gAt. 

Finally, equation ( 1153b is responsible for the change of Jo 
due to advection and compression/decompression of radiation 
together with the fluid. We "solve" this equation in the follow- 
ing way. Once the velocity-independent part of the DDMC 



20 



Abdikamalov et al. 



transport and Doppler shift operations are performed, the ad- 
vection and compression/decompression effects are modeled 
by just changing the position r of a MC particle by V,.(r) Ar„, 
where K(r) is the radial velocity of the matter at point r. This 
operation is performed at each timestep (after each DDMC 
transport step), and it automatically incorporates both advec- 
tion and compression/decompression effects. 

9. NUMERICAL IMPLEMENTATION AND TESTS 

In this Section, we describe some details of the implemen- 
tation and tests of our ID spherically symmetric code. The 
code uses spherical polar coordinates, and can handle both 
equidistant and non-equidistant grids. 

The code is currently parallelized employing hybrid 
OpenMP/MPI parallelization using the mesh replica- 
tion method (cf. Section [TJ. It uses th e open-source 
Cact u s Computational To olkit (iGoodale et al.l 
I2003L 'http://www.cactuscode.org), which provides MPI 
parallelization, input/output, and restart capability. 

9. 1. Non-moving background with fixed T and Yg 

As a first test, we consider a scenario in which radiation 
propagates through static matter with fixed temperature T and 
electron fraction Yg. 

9.1.1. Static scattering atmosphere 

iHummer & RybickH (1197 Ih found stationary-state solutions 
of the spherical analogue of the classical Milne problem. The 
model consists of a central radiation source surrounded by a 
static, spherically symmetric scattering atmosphere of some 
radius R^tm- The source emits radiation isotropically with con- 
stant luminosity Lq. The opacity of the atmosphere is assumed 
to be due only to isotropic scattering with a simple power-law 
dependence on radius: 



K, = r " , < r < /?atm , « > 1 . 



(157) 



According to the solution of iHummer & Rvbickil (Il971h . the 
luminosity L{r) at any distance from the source should be 
constant and equal to Lq. In the tests we perform, we check 
whether the condition L(r) = Lq is fulfilled for all < r < /?atm- 

We choose /?atm = 50 km and n= 1.1. The central radial 
resolution Ary is 200 m, and the entire computational domain 
is covered by A^, = 200 cells with logarithmically increasing 
size. We performed four different runs with different MCP 
weights, creating 100, 400, 1600, and 6400 new MCPs in each 
timestep for a given luminosity, Lq = 6.5 x 10'"^ erg/s 

Figure [T] shows the radial profile of the luminosity when 
a steady-state regime is reached. In agreement with the an- 
alytical solution, the luminosity indeed remains ~ Lq within 
statistical errors independent of the distance from the source. 
Such errors are expected due to the probabilistic nature of the 
Monte Carlo method. Moreover, according to the central limit 
theorem, these errors should decrease as 1 /^/N, where A^ is 
the total number of simulated MCPs. Figure |2] shows the rel- 
ative deviation AL(r)/Lo of L(r) from its correct value Lq for 
these four runs (i.e., here AL{r) = |L(r)-Lo|). As we can see, 
AL{r)/Lo indeed decreases by ~ 2 as we increase the num- 
ber of MCPs by 4. We have performed additional simulations 
with different values of «, Aro, Nr and R^tm, and the code 
again reproduces the L{r) = Lq solution with statistical errors 
that decrease as I/Vn, as in the case above. 




r [km] 

Fig. I. — The radial profile of the luminosity for a central point source 
emitting into a static scattering atmosphere for four different MC runs. In 
the first ran (black line), we emit 100 MCPs at each timestep, while in the 
second ran we produce 400 MCPs in a timestep, etc. In agreement with 
the analytical solution (cf. Section |9. Lit , the luminosity is constant along 
the radial coordinate within statistical eiTors. According to the central limit 
theorem and as demonstrated in Fig.|2] these en'ors decrease as l/N ' , where 
N is the number of MCPs. 




/[ms] 

Fig. 2. — Average deviation of the luminosity from its average value in each 
cell as a function of time for radiation from a central point source emitting 
into a static scattering atmosphere. The black line corresponds to a ran in 
which we emit 100 MCPs in each timestep, while the red line represents the 
run with 400 MCPs, etc. In agreement with the central limit theorem, the 
average deviation decreases as l/N^/~, where N is the number of MCPs. 

9.1.2. Homogeneous spiiere 

The homogeneous sphere problem is frequently em- 
ployed to test radiative transfer codes (iBruennI 119851 
ISchin der & Bludmanlll989l;ISmit et al]|1997l;lRampp & Jankal 
2002). This problem consists of a static homogeneous and 
isothermal sphere of radius R that radiates into vacuum. In- 
side the sphere, the radiation interacts with the background 
matter only via isotropic absorption and thermal emission. 

Despite its simplicity, the problem possesses some impor- 
tant physical and numerical properties that are often encoun- 
tered in many practical applications: There is a sharp dis- 
continuity at the surface of the sphere, and this represents 
a major challenge for finite-difference methods. However, 
Monte Carlo methods are well positioned for treating such 
discontinuities. Moreover, a situation similar to the stiff tran- 
sition from the radiation diffusion regime inside an opaque 
sphere to a free-streaming regime in the ambient vacuum oc- 
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curs near a PNS surface in core-collapse supernova simula- 
tions. Such a transition is a source of significant errors in 
approx imate transport schemes such as, e.g., flux-limited dif- 
fusion (lOtt et al.ll2008l) . 

We assume that the sphere of radius R has a constant ab- 
sorption opacity k,, and emissivity B in the interior, while 
in the ambient vacuum at r > R, we have Kg = B = 0. For 
this p roblem, the tran sport equation can be solved analyti- 
cally (lSmitetal.lll997h : 
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I(r,tx) = B(l-e 
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where 



r^i + Rg(r,^i) if r <R, -1 < ^t < 1 , 



s(r,ti)= I 2Rg{r,fi) if r>R, \J I - {jf < l-i < 1. 



otherwise 



and 



g(r,n)-- 



©<-"='■ 
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Note that this solution depends only on three parameters: Kg, 
R, and B, where the latter acts as a scale factor for the solution. 
We perform a set of simulations with R = 10 km. Kg = 
2.5 X 10"'* cm"' and B = 10 (in CGS units). Our computa- 
tional domain has an outer radius of 50 km and is covered 
by 100 equidistant cells. For this setup, we carry out three 
runs, in which we choose MCP weights such that 10^, 4 x 10^, 
or 1.6 X 10^ new MCPs are emitted in each timestep in a 
given simulation (hereafter, simulations A, B, and C). Fig- 
ure [3] shows the zeroth moment 7 as a function of the radial 
coordinate for simulations A, B, and C when the stationary- 
state radiation field is reached (black, red, and green lines, 
respectively). Also shown (blue line) is the zerot h moment / 
from the analytical solution given by equations ( ll58b -( fT60t . 
As we can see, the Monte Carlo solution agrees well with the 
analytical solution within statistical errors. The inset plot in 
Fig.Q]shows the time evolution of the average deviation of the 
zeroth moment in the Monte Carlo solution from the analyt- 
ical result when the stationary state is reached. The average 
deviation in simulations A, B, and C fluctuate around ~ 0.14, 
~ 0.27, ~ 0.53, respectively. Hence, in these simulations, 
when the number of MCPs is increased by a factor of 4, the 
average deviation decreases by a factor of ~ 2, in agreement 
with the central limit theorem, and the solution converges to 
the analytical result. We have repeated this simulation set with 
different values of «„, B, R and different grid resolutions. In 
all cases, we find excellent agreement with the analytical re- 
sult, similar to the case discussed above. 

9.1.3. Diffusion of a Gaussian Pulse 

In order to show that our code handles diffusion of radia- 
tion properly, we calculate the diffusion of a Gaussian radi- 
ation packet through static matter We assume that radiation 
interacts with matter only via isotropic and isoenergetic scat- 
tering, and the scattering opacity Ks is assumed to be constant 
in space and time. 

The diffusion of a Gaussian packet with initial central po- 
sition at r = rini and width do in such a medium is described 
by the following analytical solution (.Swestv & Mvra ,2009; 
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Fig. 3. — The zeroth moment for radiation / as a function of radial coor- 
dinate r in the homogeneous sphere test for three simulations in which lO'', 
4 X 10^, and 1.6 X 10* new MCPs are emitted in each timestep (black, red, 
and green lines, respectively). Also shown (blue line) is the zeroth moment 
from the analytical solution. The inset plot shows the average deviation of 
the zeroth moment in these three Monte Carlo simulations from the analyti- 
cal solution. 

ISumivoshi & Yamadall20I2l) 



E{r,t) = E, 



tini ~r t 



exp 



4D(fini + f) 



(161) 



where E{r,t) is the radiation energy density at position r and 
time t after initial time fini. The diffusion coefficient, D, equals 
c/3ks and the width of the Gaussian pulse, di^i, as a function 
of the initial time equals (4Dfini)'''^. Parameter Eini is the ini- 
tial height of the packet and uj is related to the number of 
spatial dimensions, A^d, and is equal to Nu/2. 

In our runs, we place the packet at the center of our com- 
putational domain with the radial coordinate extending up to 
3 X 10^ cm (300 km). We choose the initial width of the 
packet to be 10* cm, while the scattering opacity is set equal 
to 2 X 10""* cm"' . We run each simulation with the DDMC and 
MC schemes and use 2x10^ MCPs with constant weight in 
each of the runs. Figure |4]shows the radial profile of E at three 
different times (0, 30, and 60 ms). The black line corresponds 
to the analytical solution, the dashed-red line is obtained from 
the DDMC run, and the dotted green line is the MC solution. 
Due to diffusion, the radial profile of radiation flattens with 
time. During the first 30 ms of evolution, the central radia- 
tion energy density decreases be a factor of ^20, while in the 
next 30 ms it decreases further by a factor of ^2. Therefore, 
the 60 ms timescale captures the diffusion timescale of the 
problem. As we can see in the plot, both the DDMC and MC 
methods model diffusion of radiation in scattering medium 
quite well. The fluctuations around the analytic solution are 
due to the Monte Carlo treatment of transport in a scattering 
medium, and their magnitude again decreases as A^'/^. We 
have also repeated these runs with different values of the scat- 
tering opacity k, and different widths, din,, of the Gaussian, 
and we always find that both from DDMC and MC schemes 
agree with the analytical solution within statistical errors. 

9.2. Protoneutron star cooling 

In this section, we consider the early cooling and delep- 
tonization of a young non-rotating PNS formed in a CCSN. 
This problem is provided as a test of neutrino-matter cou- 
pling and as a testbed for optimal sampling methods and MCP 
weights. This particular problem was chosen because it is a 
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and then by choosing the weight of MCPs in each zone. We 
choose the weights of MCPs based on the following function 
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Fig. 4. — The radial profile of the energy density of a Gaussian radiation 
packet diffusing into a uniform medium with a constant scattering opacity at 
three different times. The black line coiTesponds to the analytical solution, 
the dashed red line is obtained using the DDMC scheme, while the dotted 
green line is produced by the MC scheme. Note that the fluctuations in the 
energy density that are particularly visible near the central region at ? = 30 
ms and f = 60 ms are simply caused by the random noise intrinsic to Monte 
Carlo methods. This noise is particularly pronounced in central zones be- 
cause those zones have smaller volumes and thus contain smaller numbers 
of MCPs. One can reduce this noise by choosing smaller weights for MCPs 
in central regions. In our test, we use constant weight everywhere for the 
entire radiation packet. Moreover, these fluctuation decrease with increasing 
number of MCPs, according to the central limit theorem. 

realistic physical context similar to that for which our MC and 
DDMC algorithms were designed. The reader should note, 
however, that there is no analytic or agreed-upon benchmark 
solution for this problem and, therefore, that we are testing the 
behavior and speed of the solutions, not the numerical results 
themselves. 

9.2.1. Numerical Setup 

For our PNS model, we employ the post-core bounce con- 
figuration pr oduced in the collap se of the s20.0 progeni- 
tor model of IWooslev et a J (120021) with th e 2D multi-group 
multi-angle simulations of lOtt et aP (1200811^^1 We start with a 
background model at 160 ms after bounce, and evolve it with 
our Monte Carlo code. In doing so, we first evolve the radi- 
ation field with fixed T and Y^ until it reaches a steady state, 
after which we inaugurate coupling to T and Fg- Our computa- 
tional domain extends up to 300 km, and is covered with 100 
fixed logarithmically-spaced radial zones. The central zone 
has a width of 500 m. The timestep is chosen to be the light- 
crossing time of the central zone. We assume that the PNS is 
static and neglect velocity-dependent effects. 

We include the standar d set of neutrino interactions listed 
in [Thompson et al.1 (l2003l) . We use 48 logarithmically spaced 
energy groups from 2.5 MeV to 250 MeV to calculate the 
"leakage" opacities from one energy group to another due to 
inelastic scatterings in the DDMC region (however, the ener- 
gies of MCPs are, of course, not discretized, since they are se- 
lected randomly using the local emissivity. See discussion in 
Section fTTI ). Electron neutrinos and antineutrinos are treated 
independently, while we combine heavy-lepton neutrinos (f^j, 
v^t., Vr, and Vr) together into one group. 

Sampling of Monte Carlo particles. We sample Monte 
Carlo particles in each spatial zone by calculating the num- 
ber of neutrinos emitted in each zone during each timestep 

^^ The 2D data of lOtt et al] 120081) have been angle-averaged and then 
mapped to our giid. 
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where /; is a constant parameter Terms T^mm and 7e emis are 
the rates of change of T and Y^ that would occur if there were 
only emission and no absorption. Term Kp is the mean absorp- 
tion opacity; hence, l/wpC represents a timescale for absorp- 
tion of radiation. Therefore, the quantity in expression (I162l i 
is larger in regions that emit strongly with small absorption. 
Stated differently, this quantity should have larger values in 
those regions where T and Y^ are likely to undergo significant 
changes due to emission and subsequent escape of radiation, 
implying that sampling based on quantity ( I162l i places par- 
ticular emphasize on accurate modeling of PNS cooling and 
deleptonization. We experimented with several values of h 
and found that h = 0.5 for electron neutrinos and antineutrinos 
and h = 0.3 for ji/t neutrinos lead to fairly smooth sampling 
of MCPs, where the largest number of MCPs are concentrated 
in zones that are subject to the fastest changes in T and Y^. 
However, we did not perform further studies of sampling and 
we do not claim that the spatial sampling of MCPs based on 
expression ( |162| i is the optimal choice for modeling PNS evo- 
lution. 

The Interface between DDMC and pure Monte Carlo re- 
gions. Since neutrino absorption and scattering cross sections 
are strongly energy-dependent, it is important to take this into 
account in choosing the spatial location of the interface be- 
tween the DDMC and pure MC regions. To accomplish this, 
we introduce energy groups and calculate optical depth for 
each energy group. We then determine the spatial location of 
the interface for an MCP with an energy within a given group 
in terms of the optical depth for that group. More specifically, 
if the optical depth for the energy group of a given MCP ex- 
ceeds some threshold value, tddmc, then we assume that this 
MCP is in the DDMC region. The parameter tddmc is as- 
sumed to have the same value for each energy group, mean- 
ing that the interface is located at different radii for different 
energy groups. In our simulations, we calculate the optical 
depth at each timestep and adjust the interface to its new loca- 
tion corresponding to the new values of optical depth. In our 
tests below, we consider several values of tddmc and explore 
which value represents an optimal choice for simulations of 
PNS evolution. 

Treatment of effective scattering. As mentioned in Sec- 
tion |71 in the DDMC regime, we split the effective scatter- 
ing opacity des.j.k into two parts, ajkdes.j.k and {l-ajkiaes.j.k, 
where < a^i < 1, and restrict the first of these two to be 
elastic effective scattering, while the second one is free to be 
inelastic (as "effective-scattering" scattering would otherwise 
generally be). As we have discussed in Section|7] the presence 
of an extra elastic scattering source does not increase the cost 
of doing DDMC transport. Therefore, by assuming that an 
Oj.k&es.j.k fraction of effective scattering is elastic (instead of 
being inelastic), our computational savings are proportional 
to aj,k. 

We find that it is convenient to determine the parameter aj^ 
in terms of the Fleck factor fj in the following way: 



— f I- 



(163) 



where iJ is a constant ranging from to 1 . For 3 = 0, all of ef- 
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fective scattering is treated as inelastic, while for ^ = 1, all of it 
is treated as elastic. This prescription has the advantage that at 
low optical depth - where effective scattering does not dom- 
inate calculations - most effective scattering is treated as in- 
elastic, while at high optical depth - where calculations would 
otherwise be dominated by inelastic effective scatterings - a 
significantly larger fraction is treated as elastic, which leads 
to huge savings in computation. In the following, we explore 
what values of S are most suitable for simulations of PNS evo- 
lution. 

9.2.2. Results 

In this section, we first present the stationary-state radia- 
tion field results, after which we describe the subsequent fully 
time-dependent calculations with coupling to T and Yg. For 
these simulations, unless otherwise noted, we use tddmc = 6, 
(5 = 0.38 and a = 1, and employing 100,000 MCPs to model 
newly-emitted particles at each timestep. 

Stationary state. Fi gur e |5] shows the radial profiles of the 
RMS neutrino energieofor the three types of neutrinos. The 
solid lines are produced by our Monte Carl o code, while th e 
dashed ones are from the Sn calculations of lOtt et all (120081) . 
The RMS energies agree well in the inner ^ 50 km region, 
while for r > 50 km, the Sn code produces RMS energies 
that are larger by up to ~ 4%. We believe that this difference 
st ems from the tru ncation errors in the energy discretization 
of lOtt et al] ( 120081) calculations, which employ only 16 loga- 
rithmic energy groups, whereas in our MC cases, we do not 
use energy discretization in selecting MCP energies. The cor- 
responding spectral energy distribution of the neutrino lumi- 
nosity at 300 km normalized by the total luminosity for each 
neutrino type is shown in Fig.|6] which again demonstrates a 
good a greement betwee n our MC calculation and Sn calcula- 
tions of lOtt et al.l (|2008IP I 

Figure |7] shows the radial profile of the mean inverse flux 
factoro for the three types of neutrinos. The solid lines 
again represent the Monte Carlo re sults, wh i le the dashed 
lines are from the Sn calculations of lOtt et aTl (l2008h . Over- 
all, we again find good agreement between the two results. 
The only systematic difference that we observe is that the Sn 
code tends to yield inverse flux factors that are larger by ^ 1 
% than the corresponding MC data at radii > 200 km. This is 
most likely caused by the ray-effects intrinsic to Sn schemes. 
Nevertheless, we find good overall agreement between our 
Monte Carlo and the Sn calculations of Ott et al. (2008!) for 
this stationary-state radiation field. 

Time-dependent calculations. Figure[8]shows the luminosi- 
ties for the three different neutrino species as a function of 
time during the first 150 ms of evolution of the PNS model af- 
ter being mapped to our code. As expected, the luminosities 
for all types of neutrinos decrease gradually during the first 
150 ms of evolution. On a short timescale of one timestep, 

23 See equation (13) o flOtt et al J i2008h for the exact definition of the RMS 
neutrino energy that we employ in our study. 

-■* The reason for showing normalized neutrino spectra is the following. 
We find that while the total luminosities of the heavy lepton neutrinos in 
our MC and Sn calculations agree very well, the luminosities of electron 
neutrinos and antineutrinos agree to only ~ 10-15%. Th e reason for this 
is that the angular averaging of the 2D background data of lOtt et alj J2008I) 
leads to differences in p, T, and Ye. The luminosities of electron neutrinos 
and antineutrinos are very sensitive to the values of Yc, and differences in the 
latter lead to quantitative differences in the predicted luminosities and, thus, 
to a shift in the values of the absolute spectral energy distribution. 

2^ See equation ( 1 2) of lOtt et aP J2008I) for the exact definition of the mean 
inverse flux factor that we use in our analysis. 



the luminosity undergoes significant fluctuations around its 
average value. These are expected and are due to the stochas- 
tic nature of MCP emission, absorption, and scattering pro- 
cesses. These fluctuations decrease significantly if averaged 
over longer timescales, such as the PNS light-crossing time, 
and are practically "invisible" on much longer timescales, 
such as the dynamical timescale of the PNS. Moreover, we 
find that these fluctuations again decrease as ~^ A^"'/^, where 
A^ is the number of MCPs. 

Figure |9] shows the radial profile of the temperature (upper 
panel) and electron fraction (bottom panel) at the beginning 
(t = 0) and at the end (f = 150 ms) of our simulation. During 
the early evolution, the temperature decreases noticeably in 
the region r > 15 km due to emission and diffusion of radia- 
tion. There is no significant change in T and Yf. in the inner 
r < 15 km region as the radiation diffuses out on much longer 
timescales from that region. The electron fraction decreases 
in the region r <60 km due to copious emission of electron 
neutrinos, while in the region r>60 km Fg increases as a re- 
sult of the absorption of those neutrinos. 

Quality of energy and lepton number conservation. As dis- 
cussed in Section 16.21 our IMC scheme for neutrinos does 
not make any approximations that violate energy or lepton 
number conservation. Therefore, energy and lepton number 
should be conserved up to machine precision in practical cal- 
culations using this scheme. In our calculations, we see that 
this is indeed the case, and both energy and lepton number 
in the system are conserved in each timestep to the 0(10"'^) 
precision of the machine. 

Treatment of effec tive s cattering. As we mentioned above, 
we use prescription (I163t in order to treat effective scattering 
in a computationally efficient way in the DDMC region. We 
have performed simulations for different values of the param- 
eter (5, ranging from 1 (where all effective scatterings are elas- 
tic) to 0.286 (where a significant fraction of effective scatter- 
ings are inelastic). We find that there is no systematic differ- 
ence in the PNS evolution for values of 5 that are smaller than 
~ 0.38. For example, the total energy Eesc (or lepton num- 
ber Fi.esc) of neutrinos that leave the system through the outer 
boundary in the first 150 ms of evolution decreases mono- 
tonically by ^ 2% when decreasing 5 from 1 to 0.38. Fur- 
ther decrease of 5 results in non-systematic variation in both 
/iesc and yjesc with relative error of only '--^ 0.1%, without any 
monotonic dependence on 5. The origin of such variations is 
likely to be numerical errors and not our prescription for the 
treatment of effective scatterinsPi This implies that condition 
( 1117b for the validity of our approximate treatment of effective 
scattering is fulfilled for 5 < 0.38 and violated otherwise. On 
the other hand, because the cost of simulations increases with 
decreasing 5, it is desirable to use the largest allowed value of 
5. Therefore, we conclude that 5 = 0.38 is an optimal choice 
for the treatment of effective scattering in the DDMC regime, 
at least in modeling the early phases of PNS evolution. 

Note that this feature is independent of the number of di- 
mensions. It should be applicable as long as we have a young 
PNS surrounded by a lower density envelope, i.e. in any 
CCSN scenario. 

The interface between DDMC and pure Monte Carlo re- 
gions. Figure[TO]shows the computational time to perform the 

2* Indeed, we find that these variations tend to increase when we increase 
the parameter ? that controls the re main ing fraction of an MCP when it is 
assumed to be absorbed (cf. Section |4~11 - the larger ? is. the larger are the 
numerical errors in treating MCP absoiption. 
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Fig. 5. — The stationary state radial profiles of the RMS neutrino ener- 
gies for the three different types of neutrinos obtained with our Mo nte Carlo 
schem e (solid lines). For comparison, we show results obtained by lOtt et al.l 
12008ft with an Sn code (dashed lines). 

first 500 timesteps as a function of the parameter tddmc- We 
have performed simulations using different values of tddmc, 
ranging from 5 up to 1280. Since Monte Carlo calculations 
are more expensive than DDMC calculations, the computa- 
tional time increases with tddmc- For example, for tddmc = 5 
the first 500 timesteps are performed within 1469 s, while for 
''DDMC = 1280 this simulation took 20, 480 s. We have not per- 
formed simulations with higher values of tddmc, since such 
simulations become quite expensive, implying that it is im- 
practical to perform neutrino transport in the PNS using only 
the IMC scheme without combining it with the DDMC (or 
an alternate) scheme. Interestingly, for this setup, the com- 
putational cost of simulations does not differ much for values 
of Tddmc smaller than ~ 100. This is because the DDMC 
scheme yields the largest speed-up in regimes dominated by 
scattering. At optical depth of < 100 (where the Fleck fac- 
tor / is ^ 1), there is significant contribution from absorption, 
which leads to modest increase of the computational time as 
we increase tddmc from 5 to 100. At higher optical depths 
(where / ^ 0), effective scatterings start to dominate, which 
leads to steep increases of computational time with tddmc- 

From a theoretical point of view, the DDMC scheme is ac- 
curate if the underlying approximations - which are based 
on Fick's law - are fulfilled to a sufficient degree in each 
zone. Fick's law holds in a given zone if the mean-free- 
path of MCPs are significantly smaller than the grid size. If 
we assume that the grid size is of the order of 1 km, then 
for a typical young PNS model, this requirement translates 
to Tddmc ^ 6 (for a given energy group). Therefore, for 
Tddmc ^ 6 we expect our DDMC scheme to yield sufficiently 
accurate results. Indeed, to verify this premise we have per- 
formed a series of long-evolution simulations for values of 
Tddmc in a wide range. We find that the results (such as lu- 
minosity, etc.) indeed agree to within ^ \% for any value of 
Tddmc larger than ^ 6. 

The importance of the implicit scheme. The simulations 
presented in this section are performed using the implicit 
Monte Carlo scheme with a "most implicit" choice of q = 1 
within the framework of the IMC (cf. equation l60t . We 
repeated some of the runs with the "less implicit" value of 
a = 0.5, which results in a Fleck factor that is larger by a fac- 
tor of ~ 2 than in the a = 1 case in the diffusive regime (cf. 
equation l63Tl. This results in a slight decrease in solution ac- 
curacy, but we do not observe any instabilities. One obvious 
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[in seconds] . The solid black lines show the simulation time using our code, 
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Fig . 10. — Time spent for performing first 500 timesteps as a function of pa- 
rameter tddmc for PNS evolution simulations. Since at high optical depth the 
pure Monte Carlo calculations are more expensive than the DDMC calcula- 
tions, the computational time increases with tddmc- Beyond t^dmc > 1280, 
our simulations become computationally too expensive, and we do not con- 
sider even larger values of tddmc- 



question to ask is: would it be possible to perform the same 
simulations with an explicit treatment of the emissivity? The 
answer is clearly no. With the current timestep set by the 
light crossing time for the central zone of width 500 m, the 
fully explicit scheme crashes within a few timesteps. For it to 
be stable, one would have to use timesteps of the order of the 
smallest mean-free-path, which can be as small as < 1 m in 
the center of the PNS. Hence, in order to use a time-explicit 
scheme, we would need to decrease our timestep by a factor 
of at least ^ 500, which would make useful simulations com- 
pletely impractical. For the above choice of parameters (i.e., 
TDDMC = 6, 5 = 0.38 and a = 1, and employ 100, 000 MCPs to 
model newly emitted particles at each timestep), the computa- 
tional cost of performing the first 150 ms of PNS evolution is 
about 36 hours on 96 cores on the Hopper Cray XTE6 cluster 
at NERSC. 

Parallel scaling. In order to study parallel scaling of our 
code, we perform simulations on a number of cores rang- 
ing from 24 up to 1 152 on the Hopper Cray XTE6 cluster at 
NERSC. We implement hybrid OpenMP/MPI parallelization 
and use 6 OpenMP threads per each MPI process in our tests. 

Figure [TT] shows the time spent per MPI process fnipi as 
a function of the number of cores A^core in our weak scaling 
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Fig. 12. — The radial profile of the radiation energy density measured by 
a comoving observer at different times in the homologously expanding shell 
test for simulation AID using 1,280,000 MC particles. See the discussion in 
the main text for further details of the simulation. 

tesQ In this test, we increase the number of MCPs propor- 
tionally to the number of MPI processes, while other problem 
parameters remain fixed. Here, the amount of computational 
work is chosen in such a way that the time on 24 cores is equal 
to 1 s. This time increases slowly withA^core and reaches 1.35 
s for A^core = 1 152, a value that is only 35% larger than the one 
for 24 cores. 

The inset plot of Fig. [TT]shows the total simulation time as 
a function of the number of cores in our strong scaling test. 
The simulation time decreases almost linearly (in logarithmic 
scale) with A^core and reaches 1 .53 s at A^core = 1 152. This num- 
ber is only 70% larger than the time corresponding to ideal 
scaling. Based on these results, we conclude that our code 
scales nicely up to 1152 cores, which is more than sufficient 
for ID problems. 

9.3. Moving background with fixed T and Y^ 

In this section we describe a test problem that involves the 
propagation of radiation in moving matter to demonstrate the 
ability of our code to perform transport in a moving medium. 
In this test, we neglect the change in temperature (T) and elec- 
tron fraction (Y,,) of matter due to emission, absorption, or 

^^ In the strong scaling test, the problem size remains the same as we 
increase the number of CPUs so that the amount of work per CPU decreases 
proportionately with the number of cores. Alternatively, in the weak scaling 
test, the problem size on each core remains the same, while the problem size 
increases proportionately with the number of cores used. 
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Fig . 13 . — The radiation energy density in the comoving frame as a function 
of the radial coordinate inside homologously expanding shell 

scattering of radiation. 

9.3.1. Homologously expanding shell 

A spherically-symmetric shell of matter, in which radi- 
ation is trapped, is expanding self-similarly with velocity 
Vr oc r. For simplicity, we assume that matter interacts 
with radiation only through elastic scattering, without any 
emis sion or absorption. Acco rding to the analytical solu- 
tion (iMihalas & Mihalaall984) . the average energy of neu- 
trinos (e^) in the expanding shell should decrease as 



1 



(164) 



due to the Doppler shift, while the radiation energy density £0 
should decrease according to 



Eo (X r 



(165) 



due to both the Doppler shift and expansion of the shell. 

We perform two sets of test simulations. In the first set, 
we use a very large constant value of the scattering opacity 
Ks of 5 • 10^ cm"', while in the second set, we use a million 
times smaller value of 5 • 10"^ cm"' . The reason for choosing 
two such sets is the following. Generally speaking, for such 
tests, it is desirable to have as large an opacity as possible 
in order to avoid diffusion of radiation out of the expanding 
shell, which can lead to the violation of relation ( |165t . How- 
ever, the high computational cost of performing MC trans- 
port at high scattering opacity limits the maximum value of 
Ks that we can simulate. On the other hand, the DDMC 
scheme is particularly efficient in the regime dominated by 
elastic scattering (see the discussion in Section |7]i. There- 
fore, we perform simulations with Ks = 5 x 10^ cm"' using 
only the DDMC scheme, after which we do a set of runs with 
Ks = 5x 10"-' cm"' using both the DDMC and MC schemes. 
In the former case, the radiation is completely trapped in the 
matter, and there is no diffusion of MCPs from one cell to an- 
other during the simulated time. Instead, the MCPs are trans- 
ported due only to velocity-dependent effects caused by the 
expansion of the shell. Moreover, in this case, the random 
numbers are used only for initial sampling of MCPs, while 
the velocity-dependent (operator-split) part of the transp ort i s 
completely deterministic (see the discussion in Section [8^ . 
Therefore, in this set of runs, we exclusively test the ability 
of our DDMC scheme to treat velocity-dependent effects. In 
the second set of test runs, since we use a smaller value of the 



scattering opacity, the result is affected also by diffusion of 
MCPs from the expanding shell. This is prone to statistical er- 
rors arising from the DDMC treatment of diffusion. Hence, in 
this case, we explore the ability of the code to perform trans- 
port of radiation in moving matter when the transport process 
itself is affected by statistical errors. 

In both sets of test runs, our shell initially has a radial extent 
ranging from r„iin = 2.25 km to r,nax = 12.25 km. We assume 
that the radiation initially has a constant energy density pro- 
file. The shell velocity is given by v,- = ^cr/Rout, where 03 is a 
constant, c is the speed of light, r is the radial coordinate, and 
Rout is the radius of the outer boundary of our computational 
domain. Our entire domain is covered by 1200 radial cells 
with constant size. In these simulations, we choose 03 to be 
0.4 and 0.8, while /?out is chosen to be 600 km. For the first set 
of test runs, we perform 4 simulations with 40,000, 64,000 
256,000, 1,280,000 MC particles for each of the two values 
of parameter 03. For the second set, we perform 3 simulations 
with 10,000, 40,000, and 164,000 MCPs for each value of 
03 using the MC and DDMC schemes. We use the following 
naming convention for our test simulations; the names of the 
first set of runs start with "A," while the second set starts with 
"B." The second symbol represents the value of 03 ("1" for 
03 = 0.4 and "2" for 03 = 0.8), while the last letter indicates 
whether we use the DDMC or MC scheme ("D" for DDMC 
and "M" for MC). Thus, for example, the DDMC run with 
03 = 0.4 from the first set is denoted as AID, while the similar 
run with 03 = 0.8 is called AID. 

Figure[T2]shows the radiation energy density in the comov- 
ing frame in each radial cell as a function of the radial coor- 
dinate of those cells at ms (black line), 9.2 ms (red line), 
and 19.4 ms (green line) after the start of shell expansion for 
model AID. During this evolution, the radial extent of the 
shell increases by a factor of ~ 50, which leads to a decrease 
in radiation energy density by a factor of ~5 x 10^. Due to 
homologous nature of the expansion, the radial profile of the 
energy density should stay constant. As we can see in the plot, 
our code indeed reproduces its constant radial profile within 
statistical errors. These statistical errors are caused by the ini- 
tial pseudo-random sampling of MCP positions and energies, 
and their magnitude decreases as we increase the number of 
MCPs. 

Figure [13] shows the average energy density in several out- 
ermost zones as a function of the mean radial coordinate of 
those zones for simulations AID and A2D using 1,280,000 
MCPs (red and green lines, respectively). The outermost 
zones are expanding with the fastest velocity compared to the 
other parts of the shell and, thus, are affected by the velocity- 
dependent effects to the greatest degree. This plot also shows 
energy as a function of r according to the analytical solution 
(black line). As we see, the result of these simulations agrees 
well with the analytical result during the simulated time. 

Figure [T3] also shows the energy density in several central 
zones of the expanding shell as a function of the mean radial 
coordinate of those zones for the second set of test runs. Here 
we focus on the central zones (instead of the outermost zones, 
as we did for the models of the first set) because the scattering 
opacity in these runs is relatively small. Thus, some fraction 
of the trapped radiation diffuses out of the shell during expan- 
sion. The behavior of radiation in the central zones is affected 
by this diffusion to a much smaller degree compared to that 
in the outer zones. This plot also shows the energy density 
versus the radial coordinate from the analytical solution (red 
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line). As we see, both the MC and DDMC schemes agree 
with the analytical result within the statistical errors during 
the simulated time. These statistical errors again decrease as 
we increase the number of MCPs. 

10. CONCLUSION AND FUTURE WORK 

We have generaUzed the irn p licit M onte Carlo (IMC) 
scheme of iFleck & Cummingsl (11971b an d the discrete- 
diffusi on Monte Carlo (DDMC) method of iDensmore et al] 
(120071) to energy- and velocity-dependent neutrino transport. 
While the IMC method provides a stable and accurate time 
discretization of the transport equation, the DDMC method 
increases the computational efficiency of IMC at high optical 
depths using the diffusion approximation (as is appropriate). 

In the IMC method for photons, one uses the matter energy 
equation in order to derive a time discretization of the trans- 
port equation in which the emissivity term is treated semi- 
implicitly. Since neutrino emissivity depends not only on 
energy, but also on the electron fraction, in order to derive 
the IMC equations for neutrinos with lepton number cou- 
pling, one has to use not only the energy equation, but also 
an equation for the evolution of the electron fraction Fg- Sim- 
ilar to the photon case, the IMC neutrino transport equa- 
tion has new terms that formally describe scattering of radi- 
ation. In the Monte Carlo interpretation of the IMC equa- 
tion, these scatterings model a fraction of the absorption and 
subsequent re-emission of radiation within a timestep and are 
called effective scatterings. Such scattering effectively de- 
creases the energy exchanged between matter and radiation 
within a timestep and makes the scheme unconditionally sta- 
ble. However, unlike in the photon case, the IMC scheme for 
neutrinos has two types of effective scatterings: one in which 
the total energy is conserved (as in photon IMC), and another 
for which the total lepton number is conserved during effec- 
tive scattering (this scattering type does not exist in photon 
IMC; see the discussion in Section|6|. 

In extending the gray DDMC scheme of IDensmore et al.1 
(|2007) to the energy-dependent case, we not only added dis- 
cretization in particle energy, but also proposed a new (ap- 
proximate) way of treating effective scattering in the DDMC 
regime to achieve additional speed up. We split the total ef- 
fective scattering opacity into two parts, and the first part is 
treated as inelastic (as it would otherwise be), while the other 
part is treated as elastic. This leads to extra savings in com- 
putational cost because DDMC does not need to perform any 
additional operations to model elastic scatterings (Section|7ll. 
We parametrize the fraction of the effective scattering treated 
as elastic in terms of a new variable S and show that for 
6 ~ 0.38, this approximation is excellent for modeling neu- 
trino transport in systems that consist of a newly-born PNS 
with a hot, tenuous surrounding medium. 

Furthermore, we have extended the above schemes to in- 
clude velocity-dependence. In the pure (i.e., non-DDMC) 
Mont e Carlo regime, we u se the mixed-frame form al- 
ism (iMihalas & Kleid Il982t iHubeny & Burrowsl l2007h in 
which emissivities and opacities are defined in the comoving 
frame, while transport is performed in the lab frame. This 
is easy to implement and allows one to take into account 
velocity-dependent effects correctly for arbitrarily large (rel- 
ativistic) fluid velocities, provided that the spatial resolution 



and timestep are appropriate for capturing the motion of the 
fluid in each cell. In the DDMC regime, we use a somewhat 
different approach. Since the diffusion equation is formu- 
lated most naturally in the comoving frame, we start with the 
equation for the comoving zeroth moment of radiation accu- 
rate to 0{V /c), and split it into three equations that are then 
solved in operator-split manner One part deals with diffusion 
of radiation, the second part models advection and compres- 
sion/decompression of radiation, while the third one takes into 
account the Doppler shift of radiation energy. The first equa- 
tion is identical to the DDMC equations in the zero-velocity 
case and, thus, is solved in exactly the same way, while the 
second and thi rd pa rts can be incorporated using simple oper- 
ations (Section ing . 

In order to test and validate these schemes, we have im- 
plemented them in our ID spherically-symmetric code and 
applied it to several test problems. In particular, by consid- 
ering an early evolution of a young PNS, we show that the 
IMC method allows much larger timesteps than what would 
be possible with a fully time-explicit scheme (Section 19. 2t . 
Moreover, the DDMC scheme leads to huge saving in com- 
putational cost due to the more efficient treatment of transport 
at high optical depth. Therefore, in our scheme, we use the 
DDMC method in regions with sufficiently high optical depth, 
while in the remaining part, we employ the pure IMC method. 
We find that the cost of performing a PNS evolution with the 
pure IMC scheme is usually orders of magnitude higher than 
the cost with a combination of the IMC and DDMC schemes 
(cf. Fig. ITOl i. All in all, we conclude that the combination 
of the IMC and DDMC approaches represents a robust and 
accurate method for neutrino transport in core-collapse super- 
novae. 

Finally, we point out that our energy- and velocity- 
dependent scheme can also be used for photon transport with 
minimal modifications. The only changes needed are dis- 
abling Fg changes and swapping in opacities and emissivities 
corresponding to photons. 

Having developed a velocity-dependent transport algo- 
rithm, our next task is to couple it to hydrodynamics. Fur- 
thermore, in the future, we plan to generalize our code to the 
two- and three-dimensional cases and to extend it to general 
relativity. 
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APPENDIX 
CALCULATION OF dUR/dT AND OUr/OYe USING THE FERMI-DIRAC FUNCTION 

Here, we derive analytical expressions for the equilibrium energy density of neutrinos U,- and its partial derivatives with respect 
to T and Yg using the expression for the Fermi-Dirac function B{e, T, Yg). The Fermi-Dirac function for neutrinos has the following 
form; 

B(e,r,y,)=-f-3 ^^ W^TTT' (Al) 

(hcyexp[^-r](T,Yg)\ + ^ 

where ?/ is the degeneracy parameter, which equals to (iig-ii„+^p)/kT, -{^g- ^n + jip) / kT , and for electron neutrinos, electron 
anti-neutrinos, and heavy lepton neutrinos, respectively. Here, /!«, /i„ and /i^ are the chemical potentials of the electron, neutron, 
and proton. Using the definition of the equilibrium energy density of neutrinos f/,- given by equation (|25] ), we obtain: 



dUr 


4tt 


dT 


c 


dUr 


4tt 


dYg 


c 



• dB 
dB 



BY, 



de , 



(A2) 
(A3) 



where functions |^ and ^, can be obtained using formula (lAlb : 
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Functions drj/dT and drj/dYg are expr essed in terms of the derivatives of yUp, /i„, and Hp that can be calculated using EOS tables. 
Finally, the integrals in equation (IA2b can be evaluated numerically, or can be expressed in terms of the Fermi in tegrals . The 
latter can be calculated either via direct numerical integration or by using a series expansion (iTakahashi et al.|[r978l) . 
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